diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkAddArtifactsToDwiImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkAddArtifactsToDwiImageFilter.cpp index 6fabde0dc5..21531145ef 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkAddArtifactsToDwiImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkAddArtifactsToDwiImageFilter.cpp @@ -1,337 +1,330 @@ /*=================================================================== 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 __itkAddArtifactsToDwiImageFilter_txx #define __itkAddArtifactsToDwiImageFilter_txx #include #include #include #include "itkAddArtifactsToDwiImageFilter.h" #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include namespace itk { template< class TPixelType > AddArtifactsToDwiImageFilter< TPixelType > ::AddArtifactsToDwiImageFilter() : m_UseConstantRandSeed(false) { this->SetNumberOfRequiredInputs( 1 ); m_RandGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); m_RandGen->SetSeed(); } template< class TPixelType > void AddArtifactsToDwiImageFilter< TPixelType > ::GenerateData() { if (m_UseConstantRandSeed) // always generate the same random numbers? m_RandGen->SetSeed(0); else m_RandGen->SetSeed(); m_StartTime = clock(); m_StatusText = "Starting simulation\n"; typename InputImageType::Pointer inputImage = static_cast< InputImageType* >( this->ProcessObject::GetInput(0) ); itk::ImageRegion<3> inputRegion = inputImage->GetLargestPossibleRegion(); typename itk::ImageDuplicator::Pointer duplicator = itk::ImageDuplicator::New(); duplicator->SetInputImage( inputImage ); duplicator->Update(); typename InputImageType::Pointer outputImage = duplicator->GetOutput(); // is input slize size even? int xMax=inputRegion.GetSize(0); int yMax=inputRegion.GetSize(1); if ( xMax%2 == 1 ) xMax += 1; if ( yMax%2 == 1 ) yMax += 1; // create slice object typename SliceType::Pointer slice = SliceType::New(); ImageRegion<2> sliceRegion; sliceRegion.SetSize(0, xMax); sliceRegion.SetSize(1, yMax); slice->SetLargestPossibleRegion( sliceRegion ); slice->SetBufferedRegion( sliceRegion ); slice->SetRequestedRegion( sliceRegion ); slice->Allocate(); slice->FillBuffer(0.0); ImageRegion<2> upsampledSliceRegion; if (m_Parameters.m_DoAddGibbsRinging) { upsampledSliceRegion.SetSize(0, xMax*2); upsampledSliceRegion.SetSize(1, yMax*2); } // frequency map slice - typename SliceType::Pointer fMap = NULL; + typename SliceType::Pointer fMapSlice = NULL; if (m_Parameters.m_FrequencyMap.IsNotNull()) { - fMap = SliceType::New(); - fMap->SetLargestPossibleRegion( sliceRegion ); - fMap->SetBufferedRegion( sliceRegion ); - fMap->SetRequestedRegion( sliceRegion ); - fMap->Allocate(); - fMap->FillBuffer(0.0); + fMapSlice = SliceType::New(); + fMapSlice->SetLargestPossibleRegion( sliceRegion ); + fMapSlice->SetBufferedRegion( sliceRegion ); + fMapSlice->SetRequestedRegion( sliceRegion ); + fMapSlice->Allocate(); + fMapSlice->FillBuffer(0.0); } + m_Parameters.m_SignalScale = 1; + m_Parameters.m_DoSimulateRelaxation = false; + if (m_Parameters.m_Spikes>0 || m_Parameters.m_FrequencyMap.IsNotNull() || m_Parameters.m_KspaceLineOffset>0.0 || m_Parameters.m_DoAddGibbsRinging || m_Parameters.m_EddyStrength>0 || m_Parameters.m_CroppingFactor<1.0) { ImageRegion<3> croppedRegion = inputRegion; croppedRegion.SetSize(1, croppedRegion.GetSize(1)*m_Parameters.m_CroppingFactor); itk::Point shiftedOrigin = inputImage->GetOrigin(); shiftedOrigin[1] += (inputRegion.GetSize(1)-croppedRegion.GetSize(1))*inputImage->GetSpacing()[1]/2; outputImage = InputImageType::New(); outputImage->SetSpacing( inputImage->GetSpacing() ); outputImage->SetOrigin( shiftedOrigin ); outputImage->SetDirection( inputImage->GetDirection() ); outputImage->SetLargestPossibleRegion( croppedRegion ); outputImage->SetBufferedRegion( croppedRegion ); outputImage->SetRequestedRegion( croppedRegion ); outputImage->SetVectorLength( inputImage->GetVectorLength() ); outputImage->Allocate(); typename InputImageType::PixelType temp; temp.SetSize(inputImage->GetVectorLength()); temp.Fill(0.0); outputImage->FillBuffer(temp); int tempY=croppedRegion.GetSize(1); tempY += tempY%2; croppedRegion.SetSize(1, tempY); - MatrixType transform = inputImage->GetDirection(); - for (int i=0; i<3; i++) - for (int j=0; j<3; j++) - transform[i][j] *= inputImage->GetSpacing()[j]; - m_StatusText += this->GetTime()+" > Adjusting complex signal\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"; std::vector< int > spikeVolume; for (int i=0; iGetIntegerVariate()%inputImage->GetVectorLength()); std::sort (spikeVolume.begin(), spikeVolume.end()); std::reverse (spikeVolume.begin(), spikeVolume.end()); + FiberfoxParameters doubleParam = m_Parameters.CopyParameters(); 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(inputImage->GetVectorLength()*inputRegion.GetSize(2)); for (unsigned int g=0; gGetVectorLength(); g++) { std::vector< int > spikeSlice; while (!spikeVolume.empty() && spikeVolume.back()==g) { spikeSlice.push_back(m_RandGen->GetIntegerVariate()%inputImage->GetLargestPossibleRegion().GetSize(2)); spikeVolume.pop_back(); } std::sort (spikeSlice.begin(), spikeSlice.end()); std::reverse (spikeSlice.begin(), spikeSlice.end()); for (unsigned int z=0; zGetAbortGenerateData()) { m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n"; return; } std::vector< SliceType::Pointer > compartmentSlices; // extract slice from channel g for (unsigned int y=0; yGetPixel(index3D)[g]; slice->SetPixel(index2D, pix2D); - if (fMap.IsNotNull()) - fMap->SetPixel(index2D, m_Parameters.m_FrequencyMap->GetPixel(index3D)); + if (fMapSlice.IsNotNull()) + fMapSlice->SetPixel(index2D, m_Parameters.m_FrequencyMap->GetPixel(index3D)); } if (m_Parameters.m_DoAddGibbsRinging) { itk::ResampleImageFilter::Pointer resampler = itk::ResampleImageFilter::New(); resampler->SetInput(slice); resampler->SetOutputParametersFromImage(slice); resampler->SetSize(upsampledSliceRegion.GetSize()); resampler->SetOutputSpacing(slice->GetSpacing()/2); resampler->Update(); typename SliceType::Pointer upslice = resampler->GetOutput(); compartmentSlices.push_back(upslice); - if (fMap.IsNotNull()) + if (fMapSlice.IsNotNull()) { itk::ResampleImageFilter::Pointer resampler = itk::ResampleImageFilter::New(); - resampler->SetInput(fMap); - resampler->SetOutputParametersFromImage(fMap); + resampler->SetInput(fMapSlice); + resampler->SetOutputParametersFromImage(fMapSlice); resampler->SetSize(upsampledSliceRegion.GetSize()); - resampler->SetOutputSpacing(fMap->GetSpacing()/2); + resampler->SetOutputSpacing(fMapSlice->GetSpacing()/2); resampler->Update(); - fMap = resampler->GetOutput(); + fMapSlice = resampler->GetOutput(); } } else compartmentSlices.push_back(slice); // fourier transform slice typename ComplexSliceType::Pointer fSlice; itk::Size<2> outSize; outSize.SetElement(0, xMax); outSize.SetElement(1, croppedRegion.GetSize()[1]); typename itk::KspaceImageFilter< SliceType::PixelType >::Pointer idft = itk::KspaceImageFilter< SliceType::PixelType >::New(); idft->SetUseConstantRandSeed(m_UseConstantRandSeed); + idft->SetParameters(doubleParam); idft->SetCompartmentImages(compartmentSlices); - idft->SetkOffset(m_Parameters.m_KspaceLineOffset); - idft->SettLine(m_Parameters.m_tLine); - idft->SetSimulateRelaxation(false); - idft->SetFrequencyMap(fMap); + idft->SetFrequencyMapSlice(fMapSlice); idft->SetDiffusionGradientDirection(m_Parameters.GetGradientDirection(g)); - idft->SetEddyGradientMagnitude(m_Parameters.m_EddyStrength); - idft->SetTE(m_Parameters.m_tEcho); idft->SetZ((double)z-(double)inputRegion.GetSize(2)/2.0); - idft->SetDirectionMatrix(transform); idft->SetOutSize(outSize); int numSpikes = 0; while (!spikeSlice.empty() && spikeSlice.back()==z) { numSpikes++; spikeSlice.pop_back(); } - idft->SetSpikes(numSpikes); - idft->SetSpikeAmplitude(m_Parameters.m_SpikeAmplitude); + idft->SetSpikesPerSlice(numSpikes); idft->Update(); fSlice = idft->GetOutput(); // inverse fourier transform slice typename SliceType::Pointer newSlice; typename 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++) { typename InputImageType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z; typename InputImageType::PixelType pix3D = outputImage->GetPixel(index3D); typename SliceType::IndexType index2D; index2D[0]=x; index2D[1]=y; double signal = newSlice->GetPixel(index2D); if (signal>0) signal = floor(signal+0.5); else signal = ceil(signal-0.5); pix3D[g] = signal; outputImage->SetPixel(index3D, pix3D); } ++disp; unsigned long newTick = 50*disp.count()/disp.expected_count(); for (unsigned int tick = 0; tick<(newTick-lastTick); tick++) m_StatusText += "*"; lastTick = newTick; } } m_StatusText += "\n\n"; } if (m_Parameters.m_NoiseModel!=NULL) { m_StatusText += this->GetTime()+" > Adding noise\n"; m_StatusText += "0% 10 20 30 40 50 60 70 80 90 100%\n"; m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*"; unsigned long lastTick = 0; ImageRegionIterator it1 (outputImage, outputImage->GetLargestPossibleRegion()); boost::progress_display disp(outputImage->GetLargestPossibleRegion().GetNumberOfPixels()); while(!it1.IsAtEnd()) { if (this->GetAbortGenerateData()) { m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n"; return; } ++disp; unsigned long newTick = 50*disp.count()/disp.expected_count(); for (unsigned int tick = 0; tick<(newTick-lastTick); tick++) m_StatusText += "*"; lastTick = newTick; typename InputImageType::PixelType signal = it1.Get(); m_Parameters.m_NoiseModel->AddNoise(signal); it1.Set(signal); ++it1; } m_StatusText += "\n\n"; } this->SetNthOutput(0, outputImage); m_StatusText += "Finished simulation\n"; m_StatusText += "Simulation time: "+GetTime(); } template< class TPixelType > std::string AddArtifactsToDwiImageFilter< TPixelType >::GetTime() { unsigned long total = (double)(clock() - m_StartTime)/CLOCKS_PER_SEC; 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)); return out; } } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.cpp index b1658b5849..be293d2e93 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.cpp @@ -1,223 +1,219 @@ /*=================================================================== 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 __itkKspaceImageFilter_txx #define __itkKspaceImageFilter_txx #include #include #include #include "itkKspaceImageFilter.h" #include #include #include #include #define _USE_MATH_DEFINES #include namespace itk { template< class TPixelType > KspaceImageFilter< TPixelType > ::KspaceImageFilter() - : m_tLine(1) - , m_kOffset(0) - , m_FrequencyMap(NULL) - , m_SimulateRelaxation(true) - , m_Tau(70) - , m_EddyGradientMagnitude(30) + : m_FrequencyMapSlice(NULL) , m_IsBaseline(true) - , m_SignalScale(25) - , m_Spikes(0) - , m_SpikeAmplitude(1) , m_UseConstantRandSeed(false) - , m_Tinhom(50) + , m_SpikesPerSlice(0) { m_DiffusionGradientDirection.Fill(0.0); m_RandGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); m_RandGen->SetSeed(); } template< class TPixelType > void KspaceImageFilter< TPixelType > ::BeforeThreadedGenerateData() { if (m_UseConstantRandSeed) // always generate the same random numbers? m_RandGen->SetSeed(0); else m_RandGen->SetSeed(); typename OutputImageType::Pointer outputImage = OutputImageType::New(); itk::ImageRegion<2> region; region.SetSize(0, m_OutSize[0]); region.SetSize(1, m_OutSize[1]); outputImage->SetLargestPossibleRegion( region ); outputImage->SetBufferedRegion( region ); outputImage->SetRequestedRegion( region ); outputImage->Allocate(); double gamma = 42576000; // Gyromagnetic ratio in Hz/T - if (m_EddyGradientMagnitude>0 && m_DiffusionGradientDirection.GetNorm()>0.001) + if (m_Parameters.m_EddyStrength>0 && m_DiffusionGradientDirection.GetNorm()>0.001) { - m_EddyGradientMagnitude /= 1000; // eddy gradient magnitude in T/m m_DiffusionGradientDirection.Normalize(); - m_DiffusionGradientDirection = m_DiffusionGradientDirection * m_EddyGradientMagnitude * gamma; + m_DiffusionGradientDirection = m_DiffusionGradientDirection * m_Parameters.m_EddyStrength/1000 * gamma; m_IsBaseline = false; } this->SetNthOutput(0, outputImage); m_Spike = vcl_complex(0,0); + + m_Transform = m_Parameters.m_ImageDirection; + for (int i=0; i<3; i++) + for (int j=0; j<3; j++) + m_Transform[i][j] *= m_Parameters.m_ImageSpacing[j]; } template< class TPixelType > void KspaceImageFilter< TPixelType > -::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType threadId) +::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType) { typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); typedef ImageRegionConstIterator< InputImageType > InputIteratorType; double kxMax = outputImage->GetLargestPossibleRegion().GetSize(0); // k-space size in x-direction double kyMax = outputImage->GetLargestPossibleRegion().GetSize(1); // k-space size in y-direction double xMax = m_CompartmentImages.at(0)->GetLargestPossibleRegion().GetSize(0); // scanner coverage in x-direction double yMax = m_CompartmentImages.at(0)->GetLargestPossibleRegion().GetSize(1); // scanner coverage in y-direction double numPix = kxMax*kyMax; - double dt = m_tLine/kxMax; - double fromMaxEcho = - m_tLine*kyMax/2; + double dt = m_Parameters.m_tLine/kxMax; + double fromMaxEcho = - m_Parameters.m_tLine*kyMax/2; double upsampling = xMax/kxMax; // discrepany between k-space resolution and image resolution double yMaxFov = kyMax*upsampling; // actual FOV in y-direction (in x-direction xMax==FOV) int xRingingOffset = xMax-kxMax; int yRingingOffset = yMaxFov-kyMax; while( !oit.IsAtEnd() ) { itk::Index< 2 > kIdx; kIdx[0] = oit.GetIndex()[0]; kIdx[1] = oit.GetIndex()[1]; double t = fromMaxEcho + ((double)kIdx[1]*kxMax+(double)kIdx[0])*dt; // dephasing time // rearrange slice if( kIdx[0] < kxMax/2 ) kIdx[0] = kIdx[0] + kxMax/2; else kIdx[0] = kIdx[0] - kxMax/2; if( kIdx[1] < kyMax/2 ) kIdx[1] = kIdx[1] + kyMax/2; else kIdx[1] = kIdx[1] - kyMax/2; // calculate eddy current decay factors double eddyDecay = 0; - if (m_EddyGradientMagnitude>0) - eddyDecay = exp(-(m_TE/2 + t)/m_Tau) * t/1000; + if (m_Parameters.m_EddyStrength>0) + eddyDecay = exp(-(m_Parameters.m_tEcho/2 + t)/m_Parameters.m_Tau) * t/1000; // calcualte signal relaxation factors std::vector< double > relaxFactor; - if (m_SimulateRelaxation) + if (m_Parameters.m_DoSimulateRelaxation) for (unsigned int i=0; i=kxMax/2) kx += xRingingOffset; if (ky>=kyMax/2) ky += yRingingOffset; vcl_complex s(0,0); InputIteratorType it(m_CompartmentImages.at(0), m_CompartmentImages.at(0)->GetLargestPossibleRegion() ); while( !it.IsAtEnd() ) { double x = it.GetIndex()[0]-xMax/2; double y = it.GetIndex()[1]-yMax/2; vcl_complex f(0, 0); // sum compartment signals and simulate relaxation for (unsigned int i=0; i( m_CompartmentImages.at(i)->GetPixel(it.GetIndex()) * relaxFactor.at(i) * m_SignalScale, 0); + if (m_Parameters.m_DoSimulateRelaxation) + f += std::complex( m_CompartmentImages.at(i)->GetPixel(it.GetIndex()) * relaxFactor.at(i) * m_Parameters.m_SignalScale, 0); else - f += std::complex( m_CompartmentImages.at(i)->GetPixel(it.GetIndex()) * m_SignalScale ); + f += std::complex( m_CompartmentImages.at(i)->GetPixel(it.GetIndex()) * m_Parameters.m_SignalScale ); // simulate eddy currents and other distortions double omega_t = 0; - if ( m_EddyGradientMagnitude>0 && !m_IsBaseline) + if ( m_Parameters.m_EddyStrength>0 && !m_IsBaseline) { itk::Vector< double, 3 > pos; pos[0] = x; pos[1] = y; pos[2] = m_Z; - pos = m_DirectionMatrix*pos/1000; // vector from image center to current position (in meter) + pos = m_Transform*pos/1000; // vector from image center to current position (in meter) omega_t += (m_DiffusionGradientDirection[0]*pos[0]+m_DiffusionGradientDirection[1]*pos[1]+m_DiffusionGradientDirection[2]*pos[2])*eddyDecay; } - if (m_FrequencyMap.IsNotNull()) // simulate distortions - omega_t += m_FrequencyMap->GetPixel(it.GetIndex())*t/1000; + if (m_FrequencyMapSlice.IsNotNull()) // simulate distortions + omega_t += m_FrequencyMapSlice->GetPixel(it.GetIndex())*t/1000; if (y<-yMaxFov/2) y += yMaxFov; else if (y>=yMaxFov/2) y -= yMaxFov; // actual DFT term s += f * exp( std::complex(0, 2 * M_PI * (kx*x/xMax + ky*y/yMaxFov + omega_t )) ); ++it; } s /= numPix; - if (m_Spikes>0 && sqrt(s.imag()*s.imag()+s.real()*s.real()) > sqrt(m_Spike.imag()*m_Spike.imag()+m_Spike.real()*m_Spike.real()) ) + if (m_SpikesPerSlice>0 && sqrt(s.imag()*s.imag()+s.real()*s.real()) > sqrt(m_Spike.imag()*m_Spike.imag()+m_Spike.real()*m_Spike.real()) ) m_Spike = s; outputImage->SetPixel(kIdx, s); ++oit; } } template< class TPixelType > void KspaceImageFilter< TPixelType > ::AfterThreadedGenerateData() { typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); double kxMax = outputImage->GetLargestPossibleRegion().GetSize(0); // k-space size in x-direction double kyMax = outputImage->GetLargestPossibleRegion().GetSize(1); // k-space size in y-direction - m_Spike *= m_SpikeAmplitude; + m_Spike *= m_Parameters.m_SpikeAmplitude; itk::Index< 2 > spikeIdx; - for (int i=0; iGetIntegerVariate()%(int)kxMax; spikeIdx[1] = m_RandGen->GetIntegerVariate()%(int)kyMax; outputImage->SetPixel(spikeIdx, m_Spike); } } } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.h index 63ebd14c39..cda70f2671 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.h @@ -1,134 +1,121 @@ /*=================================================================== 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 __itkKspaceImageFilter_h_ #define __itkKspaceImageFilter_h_ #include "FiberTrackingExports.h" #include #include #include #include +#include using namespace std; namespace itk{ /** * \brief Simulates k-space acquisition of one slice with a single shot EPI sequence. Enables the simulation of various effects occuring during real MR acquisitions: * - T2 signal relaxation * - Spikes * - N/2 Ghosts * - Aliasing (wrap around) * - Image distortions (off-frequency effects) * - Gibbs ringing * - Eddy current effects * Based on a discrete fourier transformation. * See "Fiberfox: Facilitating the creation of realistic white matter software phantoms" (DOI: 10.1002/mrm.25045) for details. */ template< class TPixelType > class KspaceImageFilter : public ImageSource< Image< vcl_complex< TPixelType >, 2 > > { public: typedef KspaceImageFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageSource< Image< vcl_complex< TPixelType >, 2 > > Superclass; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(KspaceImageFilter, ImageToImageFilter) typedef typename itk::Image< double, 2 > InputImageType; typedef typename InputImageType::Pointer InputImagePointerType; typedef typename Superclass::OutputImageType OutputImageType; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; typedef itk::Matrix MatrixType; typedef itk::Point Point2D; - itkSetMacro( FrequencyMap, typename InputImageType::Pointer ) ///< Used to simulate distortions. Specifies additional frequency component per voxel. - itkSetMacro( tLine, double ) ///< Time needed to fill one line in k-space (in ms). - itkSetMacro( kOffset, double ) ///< Causes N/2 ghosting artifacts. - itkSetMacro( TE, double) ///< Echo time TE (in ms). - itkSetMacro( Tinhom, double) ///< T2' signal relaxation constant (in ms). - itkSetMacro( Tau, double) ///< Eddy current decay constant (in ms). - itkSetMacro( SimulateRelaxation, bool ) ///< Enable T2 signal relaxation. + itkSetMacro( SpikesPerSlice, unsigned int ) ///< Number of spikes per slice. Corresponding parameter in fiberfox parameter object specifies the number of spikes for the whole image and can thus not be used here. + itkSetMacro( FrequencyMapSlice, typename InputImageType::Pointer ) ///< Used to simulate distortions. Specifies additional frequency component per voxel. itkSetMacro( Z, double ) ///< Slice position, necessary for eddy current simulation. - itkSetMacro( DirectionMatrix, MatrixType ) ///< Image rotation matrix, necessary for eddy current simulation. - itkSetMacro( SignalScale, double ) ///< Scaling factor for resulting signal. itkSetMacro( OutSize, itk::Size<2> ) ///< Output slice size. Can be different from input size, e.g. if Gibbs ringing is enabled. - itkSetMacro( Spikes, int ) ///< Number of randomly placed spikes per slice. - itkSetMacro( SpikeAmplitude, double ) ///< Spike amplitude relative to the largest slice value (magnitude of complex). itkSetMacro( UseConstantRandSeed, bool ) ///< Use constant seed for random generator for reproducible results. - itkSetMacro( EddyGradientMagnitude, double) ///< Magnitude of eddy current induced gradients in T/m + + void SetParameters( FiberfoxParameters param ){ m_Parameters = param; } + FiberfoxParameters GetParameters(){ return m_Parameters; } void SetCompartmentImages( std::vector< InputImagePointerType > cImgs ) { m_CompartmentImages=cImgs; } ///< One signal image per compartment. void SetT2( std::vector< double > t2Vector ) { m_T2=t2Vector; } ///< One T2 relaxation constant per compartment image. void SetDiffusionGradientDirection(itk::Vector g) { m_DiffusionGradientDirection=g; } ///< Gradient direction is needed for eddy current simulation. protected: KspaceImageFilter(); ~KspaceImageFilter() {} void BeforeThreadedGenerateData(); - void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType threadId); + void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType); void AfterThreadedGenerateData(); - bool m_SimulateRelaxation; - typename InputImageType::Pointer m_FrequencyMap; - double m_tLine; - double m_kOffset; - double m_Tinhom; - double m_TE; + FiberfoxParameters m_Parameters; + typename InputImageType::Pointer m_FrequencyMapSlice; vector< double > m_T2; vector< InputImagePointerType > m_CompartmentImages; itk::Vector m_DiffusionGradientDirection; - double m_Tau; - double m_EddyGradientMagnitude; double m_Z; - MatrixType m_DirectionMatrix; - bool m_IsBaseline; - double m_SignalScale; - itk::Size<2> m_OutSize; - int m_Spikes; - double m_SpikeAmplitude; - itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer m_RandGen; bool m_UseConstantRandSeed; + unsigned int m_SpikesPerSlice; + itk::Size<2> m_OutSize; + + bool m_IsBaseline; vcl_complex m_Spike; + MatrixType m_Transform; + itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer m_RandGen; private: }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkKspaceImageFilter.cpp" #endif #endif //__itkKspaceImageFilter_h_ diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp index 0b474f357c..f9cab62d56 100755 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp @@ -1,885 +1,867 @@ /*=================================================================== 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 namespace itk { template< class PixelType > TractsToDWIImageFilter< PixelType >::TractsToDWIImageFilter() : m_UseConstantRandSeed(false) { m_RandGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); m_RandGen->SetSeed(); } template< class PixelType > TractsToDWIImageFilter< PixelType >::~TractsToDWIImageFilter() { } template< class PixelType > TractsToDWIImageFilter< PixelType >::DoubleDwiType::Pointer TractsToDWIImageFilter< PixelType >::DoKspaceStuff( std::vector< DoubleDwiType::Pointer >& images ) { // 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(); - MatrixType transform = m_Parameters.m_ImageDirection; - for (int i=0; i<3; i++) - for (int j=0; j<3; j++) - { - if (j<2) - transform[i][j] *= m_UpsampledSpacing[j]; - else - transform[i][j] *= m_Parameters.m_ImageSpacing[j]; - } - std::vector< unsigned int > spikeVolume; for (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< 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->SetkOffset(m_Parameters.m_KspaceLineOffset); - idft->SettLine(m_Parameters.m_tLine); - idft->SetTE(m_Parameters.m_tEcho); - idft->SetTinhom(m_Parameters.m_tInhom); - idft->SetSimulateRelaxation(m_Parameters.m_DoSimulateRelaxation); - idft->SetEddyGradientMagnitude(m_Parameters.m_EddyStrength); + idft->SetParameters(m_Parameters); idft->SetZ((double)z-(double)images.at(0)->GetLargestPossibleRegion().GetSize(2)/2.0); - idft->SetDirectionMatrix(transform); - idft->SetDiffusionGradientDirection(m_Parameters.m_FiberModelList.at(0)->GetGradientDirection(g)); - idft->SetFrequencyMap(fMapSlice); - idft->SetSignalScale(m_Parameters.m_SignalScale); + 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->SetSpikes(numSpikes); - idft->SetSpikeAmplitude(m_Parameters.m_SpikeAmplitude); + idft->SetSpikesPerSlice(numSpikes); idft->Update(); ComplexSliceType::Pointer fSlice; fSlice = idft->GetOutput(); ++disp; unsigned long newTick = 50*disp.count()/disp.expected_count(); for (int 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 (int tick = 0; tick<(newTick-lastTick); tick++) m_StatusText += "*"; lastTick = newTick; } } m_StatusText += "\n\n"; return newImage; } template< class PixelType > void TractsToDWIImageFilter< PixelType >::GenerateData() { m_StartTime = clock(); m_StatusText = "Starting simulation\n"; // check input data if (m_FiberBundle.IsNull()) itkExceptionMacro("Input fiber bundle is NULL!"); if (m_Parameters.m_DoDisablePartialVolume) while (m_Parameters.m_FiberModelList.size()>1) m_Parameters.m_FiberModelList.pop_back(); 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_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); } // initialize 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); } // 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); resampler->Update(); m_Parameters.m_MaskImage = resampler->GetOutput(); } // 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); 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]; if (m_Parameters.m_DoAddMotion) { 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"; } 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 signal of " + boost::lexical_cast(m_Parameters.m_FiberModelList.size()) + " fiber compartments\n"; MITK_INFO << "Generating signal of " << m_Parameters.m_FiberModelList.size() << " fiber compartments"; int numFibers = m_FiberBundle->GetNumFibers(); boost::progress_display disp(numFibers*m_Parameters.GetNumVolumes()); ofstream logFile; logFile.open("fiberfox_motion.log"); logFile << "0 rotation: 0,0,0; translation: 0,0,0\n"; // 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(); // 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; 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(); m_StatusText += "0% 10 20 30 40 50 60 70 80 90 100%\n"; m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*"; unsigned int lastTick = 0; for (int g=0; gGetFiberPolyData(); ItkDoubleImgType::Pointer intraAxonalVolume = ItkDoubleImgType::New(); intraAxonalVolume->SetSpacing( m_UpsampledSpacing ); intraAxonalVolume->SetOrigin( m_UpsampledOrigin ); intraAxonalVolume->SetDirection( m_Parameters.m_ImageDirection ); intraAxonalVolume->SetLargestPossibleRegion( m_UpsampledImageRegion ); intraAxonalVolume->SetBufferedRegion( m_UpsampledImageRegion ); intraAxonalVolume->SetRequestedRegion( m_UpsampledImageRegion ); intraAxonalVolume->Allocate(); intraAxonalVolume->FillBuffer(0); // 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 = doubleDwi->GetPixel(idx); pix[g] += segmentVolume*m_Parameters.m_FiberModelList[k]->SimulateMeasurement(g); doubleDwi->SetPixel(idx, pix ); double vol = intraAxonalVolume->GetPixel(idx) + segmentVolume; intraAxonalVolume->SetPixel(idx, vol ); if (g==0 && vol>maxVolume) maxVolume = vol; } } ++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) fact = voxelVolume/maxVolume; while(!it3.IsAtEnd()) { if (it3.Get()>0) { DoubleDwiType::IndexType index = it3.GetIndex(); // get fiber volume fraction DoubleDwiType::Pointer fiberDwi = compartments.at(0); DoubleDwiType::PixelType fiberPix = fiberDwi->GetPixel(index); // intra axonal compartment if (fact>1) // auto scale intra-axonal if no fiber radius is specified { fiberPix[g] *= fact; fiberDwi->SetPixel(index, fiberPix); } double f = intraAxonalVolume->GetPixel(index)*fact; if (f>voxelVolume || (f>0.0 && m_Parameters.m_DoDisablePartialVolume) ) // more fiber than space in voxel? { fiberPix[g] *= voxelVolume/f; fiberDwi->SetPixel(index, fiberPix); m_VolumeFractions.at(0)->SetPixel(index, 1); } else { m_VolumeFractions.at(0)->SetPixel(index, f/voxelVolume); double nonf = voxelVolume-f; // non-fiber volume double inter = 0; if (m_Parameters.m_FiberModelList.size()>1) inter = nonf * f/voxelVolume; // inter-axonal fraction of non fiber compartment scales linearly with f double other = nonf - inter; // rest of compartment double singleinter = inter/(m_Parameters.m_FiberModelList.size()-1); // adjust non-fiber and intra-axonal signal for (unsigned int i=1; iGetPixel(index); if (f>0) pix[g] /= f; pix[g] *= singleinter; doubleDwi->SetPixel(index, pix); m_VolumeFractions.at(i)->SetPixel(index, singleinter/voxelVolume); } for (unsigned int i=0; iGetPixel(index); // if (dynamic_cast< mitk::AstroStickModel* >(m_Parameters.m_NonFiberModelList.at(i))) // { // mitk::AstroStickModel* model = dynamic_cast< mitk::AstroStickModel* >(m_Parameters.m_NonFiberModelList.at(i)); // model->SetSeed(8111984); // } pix[g] += m_Parameters.m_NonFiberModelList[i]->SimulateMeasurement(g)*other*m_Parameters.m_NonFiberModelList[i]->GetWeight(); doubleDwi->SetPixel(index, pix); m_VolumeFractions.at(i+m_Parameters.m_FiberModelList.size())->SetPixel(index, other/voxelVolume*m_Parameters.m_NonFiberModelList[i]->GetWeight()); } } } ++it3; } // move fibers if (m_Parameters.m_DoAddMotion) { if (m_Parameters.m_DoRandomizeMotion) { fiberBundleTransformed = fiberBundle->GetDeepCopy(); 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 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]); } } 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_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) { 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; } else { m_StatusText += this->GetTime()+" > Summing compartments\n"; MITK_INFO << "Summing compartments"; doubleOutImage = compartments.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->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 (int 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; } 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) 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(); } 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 > std::string TractsToDWIImageFilter< PixelType >::GetTime() { unsigned long total = (double)(clock() - m_StartTime)/CLOCKS_PER_SEC; 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)); return out; } } diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.cpp index 3314f2fc3c..4069018771 100644 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.cpp +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.cpp @@ -1,450 +1,451 @@ /*=================================================================== 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 #include #include #include "mitkFiberfoxParameters.h" #include #include #include #include #include template< class ScalarType > mitk::FiberfoxParameters< ScalarType >::FiberfoxParameters() : m_DoAddGibbsRinging(false) , m_ArtifactModelString("") , m_AxonRadius(0) , m_Bvalue(1000) , m_DoAddMotion(false) , m_DoDisablePartialVolume(false) , m_DoSimulateRelaxation(true) , m_EddyStrength(0) + , m_Tau(70) , m_KspaceLineOffset(0) , m_NumGradients(6) , m_NumBaseline(1) , m_OutputPath("") , m_DoRandomizeMotion(true) , m_Repetitions(1) , m_SignalModelString("") , m_SignalScale(100) , m_SpikeAmplitude(1) , m_Spikes(0) , m_tEcho(100) , m_tInhom(50) , m_tLine(1) , m_CroppingFactor(1) , m_MaskImage(NULL) , m_FrequencyMap(NULL) , m_NoiseModel(NULL) { m_ImageDirection.SetIdentity(); m_ImageOrigin.Fill(0.0); m_ImageRegion.SetSize(0, 11); m_ImageRegion.SetSize(1, 11); m_ImageRegion.SetSize(2, 3); m_ImageSpacing.Fill(2.0); m_Translation.Fill(0.0); m_Rotation.Fill(0.0); m_ResultNode = mitk::DataNode::New(); m_ParentNode = NULL; GenerateGradientHalfShell(); } template< class ScalarType > mitk::FiberfoxParameters< ScalarType >::~FiberfoxParameters() { // if (m_NoiseModel!=NULL) // delete m_NoiseModel; } template< class ScalarType > void mitk::FiberfoxParameters< ScalarType >::GenerateGradientHalfShell() { int NPoints = 2*m_NumGradients; m_GradientDirections.clear(); m_NumBaseline = NPoints/20; if (m_NumBaseline==0) m_NumBaseline=1; GradientType g; g.Fill(0.0); for (unsigned int i=0; i theta; theta.set_size(NPoints); vnl_vector phi; phi.set_size(NPoints); double C = sqrt(4*M_PI); phi(0) = 0.0; phi(NPoints-1) = 0.0; for(int i=0; i0 && i std::vector< int > mitk::FiberfoxParameters< ScalarType >::GetBaselineIndices() { std::vector< int > result; for( unsigned int i=0; im_GradientDirections.size(); i++) if (m_GradientDirections.at(i).GetNorm()<0.0001) result.push_back(i); return result; } template< class ScalarType > unsigned int mitk::FiberfoxParameters< ScalarType >::GetFirstBaselineIndex() { for( unsigned int i=0; im_GradientDirections.size(); i++) if (m_GradientDirections.at(i).GetNorm()<0.0001) return i; return -1; } template< class ScalarType > bool mitk::FiberfoxParameters< ScalarType >::IsBaselineIndex(unsigned int idx) { if (m_GradientDirections.size()>idx && m_GradientDirections.at(idx).GetNorm()<0.0001) return true; return false; } template< class ScalarType > unsigned int mitk::FiberfoxParameters< ScalarType >::GetNumWeightedVolumes() { return m_NumGradients; } template< class ScalarType > unsigned int mitk::FiberfoxParameters< ScalarType >::GetNumBaselineVolumes() { return m_NumBaseline; } template< class ScalarType > unsigned int mitk::FiberfoxParameters< ScalarType >::GetNumVolumes() { return m_GradientDirections.size(); } template< class ScalarType > typename mitk::FiberfoxParameters< ScalarType >::GradientListType mitk::FiberfoxParameters< ScalarType >::GetGradientDirections() { return m_GradientDirections; } template< class ScalarType > typename mitk::FiberfoxParameters< ScalarType >::GradientType mitk::FiberfoxParameters< ScalarType >::GetGradientDirection(unsigned int i) { if (i void mitk::FiberfoxParameters< ScalarType >::SetNumWeightedGradients(int numGradients) { m_NumGradients = numGradients; GenerateGradientHalfShell(); } template< class ScalarType > void mitk::FiberfoxParameters< ScalarType >::SetGradienDirections(GradientListType gradientList) { m_GradientDirections = gradientList; m_NumGradients = 0; m_NumBaseline = 0; for( unsigned int i=0; im_GradientDirections.size(); i++) { if (m_GradientDirections.at(i).GetNorm()>0.0001) m_NumGradients++; else m_NumBaseline++; } } template< class ScalarType > void mitk::FiberfoxParameters< ScalarType >::SetGradienDirections(mitk::DiffusionImage::GradientDirectionContainerType::Pointer gradientList) { m_NumGradients = 0; m_NumBaseline = 0; m_GradientDirections.clear(); for( unsigned int i=0; iSize(); i++) { GradientType g; g[0] = gradientList->at(i)[0]; g[1] = gradientList->at(i)[1]; g[2] = gradientList->at(i)[2]; m_GradientDirections.push_back(g); if (m_GradientDirections.at(i).GetNorm()>0.0001) m_NumGradients++; else m_NumBaseline++; } } template< class ScalarType > void mitk::FiberfoxParameters< ScalarType >::LoadParameters(string filename) { boost::property_tree::ptree parameters; boost::property_tree::xml_parser::read_xml(filename, parameters); m_FiberModelList.clear(); m_NonFiberModelList.clear(); if (m_NoiseModel!=NULL) delete m_NoiseModel; BOOST_FOREACH( boost::property_tree::ptree::value_type const& v1, parameters.get_child("fiberfox") ) { if( v1.first == "image" ) { m_ImageRegion.SetSize(0, v1.second.get("basic.size.x")); m_ImageRegion.SetSize(1, v1.second.get("basic.size.y")); m_ImageRegion.SetSize(2, v1.second.get("basic.size.z")); m_ImageSpacing[0] = v1.second.get("basic.spacing.x"); m_ImageSpacing[1] = v1.second.get("basic.spacing.y"); m_ImageSpacing[2] = v1.second.get("basic.spacing.z"); m_NumGradients = v1.second.get("basic.numgradients"); GenerateGradientHalfShell(); m_Bvalue = v1.second.get("basic.bvalue"); m_Repetitions = v1.second.get("repetitions"); m_SignalScale = v1.second.get("signalScale"); m_tEcho = v1.second.get("tEcho"); m_tLine = v1.second.get("tLine"); m_tInhom = v1.second.get("tInhom"); m_AxonRadius = v1.second.get("axonRadius"); m_DoSimulateRelaxation = v1.second.get("doSimulateRelaxation"); m_DoDisablePartialVolume = v1.second.get("doDisablePartialVolume"); if (v1.second.get("artifacts.addnoise")) { switch (v1.second.get("artifacts.noisedistribution")) { case 0: m_NoiseModel = new mitk::RicianNoiseModel< ScalarType >(); break; case 1: m_NoiseModel = new mitk::ChiSquareNoiseModel< ScalarType >(); break; default: m_NoiseModel = new mitk::RicianNoiseModel< ScalarType >(); } m_NoiseModel->SetNoiseVariance(v1.second.get("artifacts.noisevariance")); } m_KspaceLineOffset = v1.second.get("artifacts.m_KspaceLineOffset"); m_CroppingFactor = (100-v1.second.get("artifacts.aliasingfactor"))/100; m_Spikes = v1.second.get("artifacts.m_Spikesnum"); m_SpikeAmplitude = v1.second.get("artifacts.m_Spikesscale"); m_EddyStrength = v1.second.get("artifacts.m_EddyStrength"); m_DoAddGibbsRinging = v1.second.get("artifacts.addringing"); m_DoAddMotion = v1.second.get("artifacts.m_DoAddMotion"); m_DoRandomizeMotion = v1.second.get("artifacts.m_RandomMotion"); m_Translation[0] = v1.second.get("artifacts.m_Translation0"); m_Translation[1] = v1.second.get("artifacts.m_Translation1"); m_Translation[2] = v1.second.get("artifacts.m_Translation2"); m_Rotation[0] = v1.second.get("artifacts.m_Rotation0"); m_Rotation[1] = v1.second.get("artifacts.m_Rotation1"); m_Rotation[2] = v1.second.get("artifacts.m_Rotation2"); // compartment 1 switch (v1.second.get("compartment1.index")) { case 0: mitk::StickModel* stickModel = new mitk::StickModel(); stickModel->SetGradientList(m_GradientDirections); stickModel->SetBvalue(m_Bvalue); stickModel->SetDiffusivity(v1.second.get("compartment1.stick.d")); stickModel->SetT2(v1.second.get("compartment1.stick.t2")); m_FiberModelList.push_back(stickModel); break; case 1: mitk::TensorModel* zeppelinModel = new mitk::TensorModel(); zeppelinModel->SetGradientList(m_GradientDirections); zeppelinModel->SetBvalue(m_Bvalue); zeppelinModel->SetDiffusivity1(v1.second.get("compartment1.zeppelin.d1")); zeppelinModel->SetDiffusivity2(v1.second.get("compartment1.zeppelin.d2")); zeppelinModel->SetDiffusivity3(v1.second.get("compartment1.zeppelin.d2")); zeppelinModel->SetT2(v1.second.get("compartment1.zeppelin.t2")); m_FiberModelList.push_back(zeppelinModel); break; case 2: mitk::TensorModel* tensorModel = new mitk::TensorModel(); tensorModel->SetGradientList(m_GradientDirections); tensorModel->SetBvalue(m_Bvalue); tensorModel->SetDiffusivity1(v1.second.get("compartment1.tensor.d1")); tensorModel->SetDiffusivity2(v1.second.get("compartment1.tensor.d2")); tensorModel->SetDiffusivity3(v1.second.get("compartment1.tensor.d3")); tensorModel->SetT2(v1.second.get("compartment1.tensor.t2")); m_FiberModelList.push_back(tensorModel); break; } // compartment 2 switch (v1.second.get("compartment2.index")) { case 0: mitk::StickModel* stickModel = new mitk::StickModel(); stickModel->SetGradientList(m_GradientDirections); stickModel->SetBvalue(m_Bvalue); stickModel->SetDiffusivity(v1.second.get("compartment2.stick.d")); stickModel->SetT2(v1.second.get("compartment2.stick.t2")); m_FiberModelList.push_back(stickModel); break; case 1: mitk::TensorModel* zeppelinModel = new mitk::TensorModel(); zeppelinModel->SetGradientList(m_GradientDirections); zeppelinModel->SetBvalue(m_Bvalue); zeppelinModel->SetDiffusivity1(v1.second.get("compartment2.zeppelin.d1")); zeppelinModel->SetDiffusivity2(v1.second.get("compartment2.zeppelin.d2")); zeppelinModel->SetDiffusivity3(v1.second.get("compartment2.zeppelin.d2")); zeppelinModel->SetT2(v1.second.get("compartment2.zeppelin.t2")); m_FiberModelList.push_back(zeppelinModel); break; case 2: mitk::TensorModel* tensorModel = new mitk::TensorModel(); tensorModel->SetGradientList(m_GradientDirections); tensorModel->SetBvalue(m_Bvalue); tensorModel->SetDiffusivity1(v1.second.get("compartment2.tensor.d1")); tensorModel->SetDiffusivity2(v1.second.get("compartment2.tensor.d2")); tensorModel->SetDiffusivity3(v1.second.get("compartment2.tensor.d3")); tensorModel->SetT2(v1.second.get("compartment2.tensor.t2")); m_FiberModelList.push_back(tensorModel); break; } // compartment 3 switch (v1.second.get("compartment3.index")) { case 0: mitk::BallModel* ballModel = new mitk::BallModel(); ballModel->SetGradientList(m_GradientDirections); ballModel->SetBvalue(m_Bvalue); ballModel->SetDiffusivity(v1.second.get("compartment3.ball.d")); ballModel->SetT2(v1.second.get("compartment3.ball.t2")); ballModel->SetWeight(v1.second.get("compartment3.weight")); m_NonFiberModelList.push_back(ballModel); break; case 1: mitk::AstroStickModel* astrosticksModel = new mitk::AstroStickModel(); astrosticksModel->SetGradientList(m_GradientDirections); astrosticksModel->SetBvalue(m_Bvalue); astrosticksModel->SetDiffusivity(v1.second.get("compartment3.astrosticks.d")); astrosticksModel->SetT2(v1.second.get("compartment3.astrosticks.t2")); astrosticksModel->SetRandomizeSticks(v1.second.get("compartment3.astrosticks.randomize")); astrosticksModel->SetWeight(v1.second.get("compartment3.weight")); m_NonFiberModelList.push_back(astrosticksModel); break; case 2: mitk::DotModel* dotModel = new mitk::DotModel(); dotModel->SetGradientList(m_GradientDirections); dotModel->SetT2(v1.second.get("compartment3.dot.t2")); dotModel->SetWeight(v1.second.get("compartment3.weight")); m_NonFiberModelList.push_back(dotModel); break; } // compartment 4 switch (v1.second.get("compartment4.index")) { case 0: mitk::BallModel* ballModel = new mitk::BallModel(); ballModel->SetGradientList(m_GradientDirections); ballModel->SetBvalue(m_Bvalue); ballModel->SetDiffusivity(v1.second.get("compartment4.ball.d")); ballModel->SetT2(v1.second.get("compartment4.ball.t2")); ballModel->SetWeight(v1.second.get("compartment4.weight")); m_NonFiberModelList.push_back(ballModel); break; case 1: mitk::AstroStickModel* astrosticksModel = new mitk::AstroStickModel(); astrosticksModel->SetGradientList(m_GradientDirections); astrosticksModel->SetBvalue(m_Bvalue); astrosticksModel->SetDiffusivity(v1.second.get("compartment4.astrosticks.d")); astrosticksModel->SetT2(v1.second.get("compartment4.astrosticks.t2")); astrosticksModel->SetRandomizeSticks(v1.second.get("compartment4.astrosticks.randomize")); astrosticksModel->SetWeight(v1.second.get("compartment4.weight")); m_NonFiberModelList.push_back(astrosticksModel); break; case 2: mitk::DotModel* dotModel = new mitk::DotModel(); dotModel->SetGradientList(m_GradientDirections); dotModel->SetT2(v1.second.get("compartment4.dot.t2")); dotModel->SetWeight(v1.second.get("compartment4.weight")); m_NonFiberModelList.push_back(dotModel); break; } } } } template< class ScalarType > void mitk::FiberfoxParameters< ScalarType >::PrintSelf() { MITK_INFO << "m_ImageRegion: " << m_ImageRegion; MITK_INFO << "m_ImageSpacing: " << m_ImageSpacing; MITK_INFO << "m_ImageOrigin: " << m_ImageOrigin; MITK_INFO << "m_ImageDirection: " << m_ImageDirection; MITK_INFO << "m_NumGradients: " << m_NumGradients; MITK_INFO << "m_Bvalue: " << m_Bvalue; MITK_INFO << "m_Repetitions: " << m_Repetitions; MITK_INFO << "m_SignalScale: " << m_SignalScale; MITK_INFO << "m_tEcho: " << m_tEcho; MITK_INFO << "m_tLine: " << m_tLine; MITK_INFO << "m_tInhom: " << m_tInhom; MITK_INFO << "m_AxonRadius: " << m_AxonRadius; MITK_INFO << "m_KspaceLineOffset: " << m_KspaceLineOffset; MITK_INFO << "m_AddGibbsRinging: " << m_DoAddGibbsRinging; MITK_INFO << "m_EddyStrength: " << m_EddyStrength; MITK_INFO << "m_Spikes: " << m_Spikes; MITK_INFO << "m_SpikeAmplitude: " << m_SpikeAmplitude; - MITK_INFO << "m_Wrap: " << m_CroppingFactor; + MITK_INFO << "m_CroppingFactor: " << m_CroppingFactor; MITK_INFO << "m_DoSimulateRelaxation: " << m_DoSimulateRelaxation; MITK_INFO << "m_DoDisablePartialVolume: " << m_DoDisablePartialVolume; MITK_INFO << "m_DoAddMotion: " << m_DoAddMotion; MITK_INFO << "m_RandomMotion: " << m_DoRandomizeMotion; MITK_INFO << "m_Translation: " << m_Translation; MITK_INFO << "m_Rotation: " << m_Rotation; MITK_INFO << "m_SignalModelString: " << m_SignalModelString; MITK_INFO << "m_ArtifactModelString: " << m_ArtifactModelString; MITK_INFO << "m_OutputPath: " << m_OutputPath; } diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h index 6986f45286..d6a2e1f658 100644 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h @@ -1,175 +1,176 @@ /*=================================================================== 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 _MITK_FiberfoxParameters_H #define _MITK_FiberfoxParameters_H #include #include #include #include #include #include #include #include using namespace std; namespace mitk { /** * \brief Datastructure to manage the Fiberfox signal generation parameters. * */ template< class ScalarType > class FiberfoxParameters { public: typedef itk::Image ItkDoubleImgType; typedef itk::Image ItkUcharImgType; typedef std::vector< DiffusionSignalModel* > DiffusionModelListType; typedef DiffusionSignalModel::GradientListType GradientListType; typedef DiffusionSignalModel::GradientType GradientType; typedef DiffusionNoiseModel NoiseModelType; typedef DiffusionSignalModel* DiffusionModelType; FiberfoxParameters(); ~FiberfoxParameters(); /** Get same parameter object with different template parameter */ template< class OutType > FiberfoxParameters< OutType > CopyParameters() { FiberfoxParameters< OutType > out; out.m_ImageRegion = m_ImageRegion; out.m_ImageSpacing = m_ImageSpacing; out.m_ImageOrigin = m_ImageOrigin; out.m_ImageDirection = m_ImageDirection; out.SetNumWeightedGradients(m_NumGradients); out.m_Bvalue = m_Bvalue; out.m_Repetitions = m_Repetitions; out.m_SignalScale = m_SignalScale; out.m_tEcho = m_tEcho; out.m_tLine = m_tLine; out.m_tInhom = m_tInhom; out.m_AxonRadius = m_AxonRadius; out.m_KspaceLineOffset = m_KspaceLineOffset; out.m_DoAddGibbsRinging = m_DoAddGibbsRinging; out.m_EddyStrength = m_EddyStrength; out.m_Spikes = m_Spikes; out.m_SpikeAmplitude = m_SpikeAmplitude; out.m_CroppingFactor = m_CroppingFactor; out.m_DoSimulateRelaxation = m_DoSimulateRelaxation; out.m_DoDisablePartialVolume = m_DoDisablePartialVolume; out.m_DoAddMotion = m_DoAddMotion; out.m_DoRandomizeMotion = m_DoRandomizeMotion; out.m_Translation = m_Translation; out.m_Rotation = m_Rotation; if (m_NoiseModel!=NULL) { if (dynamic_cast*>(m_NoiseModel)) out.m_NoiseModel = new mitk::RicianNoiseModel(); else if (dynamic_cast*>(m_NoiseModel)) out.m_NoiseModel = new mitk::ChiSquareNoiseModel(); out.m_NoiseModel->SetNoiseVariance(m_NoiseModel->GetNoiseVariance()); } out.m_FrequencyMap = m_FrequencyMap; out.m_MaskImage = m_MaskImage; out.m_ResultNode = m_ResultNode; out.m_ParentNode = m_ParentNode; out.m_SignalModelString = m_SignalModelString; out.m_ArtifactModelString = m_ArtifactModelString; out.m_OutputPath = m_OutputPath; return out; } /** Output image specifications */ itk::ImageRegion<3> m_ImageRegion; ///< Image size. itk::Vector m_ImageSpacing; ///< Image voxel size. itk::Point m_ImageOrigin; ///< Image origin. itk::Matrix m_ImageDirection; ///< Image rotation matrix. /** Other acquisitions parameters */ unsigned int m_Repetitions; ///< Noise will be summed N times and afterwards averaged. double m_SignalScale; ///< Scaling factor for output signal (before noise is added). double m_tEcho; ///< Echo time TE. double m_tLine; ///< k-space line readout time. double m_tInhom; ///< T2' double m_Bvalue; /** Signal generation */ DiffusionModelListType m_FiberModelList; ///< Intra- and inter-axonal compartments. DiffusionModelListType m_NonFiberModelList; ///< Extra-axonal compartments. double m_AxonRadius; ///< Determines compartment volume fractions (0 == automatic axon radius estimation) /** Artifacts */ int m_Spikes; ///< Number of spikes randomly appearing in the image double m_SpikeAmplitude; ///< amplitude of spikes relative to the largest signal intensity (magnitude of complex) double m_KspaceLineOffset; ///< Causes N/2 ghosts. Larger offset means stronger ghost. - double m_EddyStrength; ///< Strength of eddy current induced gradients in T/m. + double m_EddyStrength; ///< Strength of eddy current induced gradients in mT/m. + double m_Tau; ///< Eddy current decay constant (in ms) double m_CroppingFactor; ///< FOV size in y-direction is multiplied by this factor. Causes aliasing artifacts. bool m_DoAddGibbsRinging; ///< Add Gibbs ringing artifact bool m_DoSimulateRelaxation; ///< Add T2 relaxation effects bool m_DoDisablePartialVolume; ///< Disable partial volume effects. Each voxel is either all fiber or all non-fiber. bool m_DoAddMotion; ///< Enable motion artifacts. bool m_DoRandomizeMotion; ///< Toggles between random and linear motion. itk::Vector m_Translation; ///< Maximum translational motion. itk::Vector m_Rotation; ///< Maximum rotational motion. NoiseModelType* m_NoiseModel; ///< If != NULL, noise is added to the image. ItkDoubleImgType::Pointer m_FrequencyMap; ///< If != NULL, distortions are added to the image using this frequency map. ItkUcharImgType::Pointer m_MaskImage; ///< Signal is only genrated inside of the mask image. /** Output parameters (only relevant in GUI application) */ mitk::DataNode::Pointer m_ResultNode; ///< Stores resulting image. mitk::DataNode::Pointer m_ParentNode; ///< Parent node or result node. string m_SignalModelString; ///< Appendet to the name of the result node string m_ArtifactModelString; ///< Appendet to the name of the result node string m_OutputPath; ///< Image is automatically saved to the specified folder after simulation is finished. void PrintSelf(); ///< Print parameters to stdout. void LoadParameters(string filename); ///< Load image generation parameters from .ffp file. void GenerateGradientHalfShell(); ///< Generates half shell of gradient directions (with m_NumGradients non-zero directions) std::vector< int > GetBaselineIndices(); unsigned int GetFirstBaselineIndex(); bool IsBaselineIndex(unsigned int idx); unsigned int GetNumWeightedVolumes(); unsigned int GetNumBaselineVolumes(); unsigned int GetNumVolumes(); GradientListType GetGradientDirections(); GradientType GetGradientDirection(unsigned int i); void SetNumWeightedGradients(int numGradients); ///< Automaticall calls GenerateGradientHalfShell() afterwards. void SetGradienDirections(GradientListType gradientList); void SetGradienDirections(mitk::DiffusionImage::GradientDirectionContainerType::Pointer gradientList); protected: unsigned int m_NumBaseline; ///< Number of non-diffusion-weighted image volumes. unsigned int m_NumGradients; ///< Number of diffusion-weighted image volumes. GradientListType m_GradientDirections; ///< Total number of image volumes. }; } #include "mitkFiberfoxParameters.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxAddArtifactsToDwiTest.cpp b/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxAddArtifactsToDwiTest.cpp index 44a19dab9f..550efe2c27 100644 --- a/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxAddArtifactsToDwiTest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxAddArtifactsToDwiTest.cpp @@ -1,193 +1,205 @@ /*=================================================================== 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 #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include /**Documentation * Test the Fiberfox simulation functions (diffusion weighted image -> diffusion weighted image) */ class mitkFiberfoxAddArtifactsToDwiTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkFiberfoxAddArtifactsToDwiTestSuite); MITK_TEST(Spikes); MITK_TEST(GibbsRinging); MITK_TEST(Ghost); MITK_TEST(Aliasing); MITK_TEST(Eddy); MITK_TEST(RicianNoise); MITK_TEST(ChiSquareNoise); MITK_TEST(Distortions); CPPUNIT_TEST_SUITE_END(); private: mitk::DiffusionImage::Pointer m_InputDwi; FiberfoxParameters m_Parameters; public: void setUp() { RegisterDiffusionCoreObjectFactory(); // reference files m_InputDwi = dynamic_cast*>(mitk::IOUtil::LoadDataNode(GetTestDataFilePath("DiffusionImaging/Fiberfox/StickBall_RELAX.dwi"))->GetData()); // parameter setup m_Parameters = FiberfoxParameters(); m_Parameters.m_ImageRegion = m_InputDwi->GetVectorImage()->GetLargestPossibleRegion(); m_Parameters.m_ImageSpacing = m_InputDwi->GetVectorImage()->GetSpacing(); m_Parameters.m_ImageOrigin = m_InputDwi->GetVectorImage()->GetOrigin(); m_Parameters.m_ImageDirection = m_InputDwi->GetVectorImage()->GetDirection(); m_Parameters.m_Bvalue = m_InputDwi->GetB_Value(); m_Parameters.SetGradienDirections(m_InputDwi->GetDirections()); } bool CompareDwi(itk::VectorImage< short, 3 >* dwi1, itk::VectorImage< short, 3 >* dwi2) { typedef itk::VectorImage< short, 3 > DwiImageType; try{ itk::ImageRegionIterator< DwiImageType > it1(dwi1, dwi1->GetLargestPossibleRegion()); itk::ImageRegionIterator< DwiImageType > it2(dwi2, dwi2->GetLargestPossibleRegion()); while(!it1.IsAtEnd()) { if (it1.Get()!=it2.Get()) return false; ++it1; ++it2; } } catch(...) { return false; } return true; } void StartSimulation(string testFileName) { mitk::DiffusionImage::Pointer refImage = NULL; if (!testFileName.empty()) CPPUNIT_ASSERT(refImage = dynamic_cast*>(mitk::IOUtil::LoadDataNode(testFileName)->GetData())); itk::AddArtifactsToDwiImageFilter< short >::Pointer artifactsToDwiFilter = itk::AddArtifactsToDwiImageFilter< short >::New(); artifactsToDwiFilter->SetUseConstantRandSeed(true); artifactsToDwiFilter->SetInput(m_InputDwi->GetVectorImage()); artifactsToDwiFilter->SetParameters(m_Parameters); CPPUNIT_ASSERT_NO_THROW(artifactsToDwiFilter->Update()); mitk::DiffusionImage::Pointer testImage = mitk::DiffusionImage::New(); testImage->SetVectorImage( artifactsToDwiFilter->GetOutput() ); testImage->SetB_Value(m_Parameters.m_Bvalue); testImage->SetDirections(m_Parameters.GetGradientDirections()); testImage->InitializeFromVectorImage(); if (refImage.IsNotNull()) { - CPPUNIT_ASSERT_MESSAGE(testFileName, CompareDwi(testImage->GetVectorImage(), refImage->GetVectorImage())); + bool ok = CompareDwi(testImage->GetVectorImage(), refImage->GetVectorImage()); + if (!ok) + { + NrrdDiffusionImageWriter::Pointer writer = NrrdDiffusionImageWriter::New(); + writer->SetFileName("/tmp/test2.dwi"); + writer->SetInput(testImage); + writer->Update(); + + writer->SetFileName("/tmp/ref2.dwi"); + writer->SetInput(refImage); + writer->Update(); + } + CPPUNIT_ASSERT_MESSAGE(testFileName, ok); } else { NrrdDiffusionImageWriter::Pointer writer = NrrdDiffusionImageWriter::New(); writer->SetFileName("/local/distortions2.dwi"); writer->SetInput(testImage); writer->Update(); } } void Spikes() { m_Parameters.m_Spikes = 5; m_Parameters.m_SpikeAmplitude = 1; StartSimulation( GetTestDataFilePath("DiffusionImaging/Fiberfox/spikes2.dwi") ); } void GibbsRinging() { m_Parameters.m_DoAddGibbsRinging = true; StartSimulation( GetTestDataFilePath("DiffusionImaging/Fiberfox/gibbsringing2.dwi") ); } void Ghost() { m_Parameters.m_KspaceLineOffset = 0.25; StartSimulation( GetTestDataFilePath("DiffusionImaging/Fiberfox/ghost2.dwi") ); } void Aliasing() { m_Parameters.m_CroppingFactor = 0.4; StartSimulation( GetTestDataFilePath("DiffusionImaging/Fiberfox/aliasing2.dwi") ); } void Eddy() { m_Parameters.m_EddyStrength = 0.05; StartSimulation( GetTestDataFilePath("DiffusionImaging/Fiberfox/eddy2.dwi") ); } void RicianNoise() { mitk::RicianNoiseModel* ricianNoiseModel = new mitk::RicianNoiseModel(); ricianNoiseModel->SetNoiseVariance(1000000); ricianNoiseModel->SetSeed(0); m_Parameters.m_NoiseModel = ricianNoiseModel; // StartSimulation( GetTestDataFilePath("DiffusionImaging/Fiberfox/riciannoise2.dwi") ); delete m_Parameters.m_NoiseModel; } void ChiSquareNoise() { mitk::ChiSquareNoiseModel* chiSquareNoiseModel = new mitk::ChiSquareNoiseModel(); chiSquareNoiseModel->SetNoiseVariance(1000000); chiSquareNoiseModel->SetSeed(0); m_Parameters.m_NoiseModel = chiSquareNoiseModel; // StartSimulation( GetTestDataFilePath("DiffusionImaging/Fiberfox/chisquarenoise2.dwi") ); delete m_Parameters.m_NoiseModel; } void Distortions() { mitk::Image::Pointer mitkFMap = dynamic_cast(mitk::IOUtil::LoadDataNode( GetTestDataFilePath("DiffusionImaging/Fiberfox/Fieldmap.nrrd") )->GetData()); typedef itk::Image ItkDoubleImgType; ItkDoubleImgType::Pointer fMap = ItkDoubleImgType::New(); mitk::CastToItkImage(mitkFMap, fMap); m_Parameters.m_FrequencyMap = fMap; StartSimulation( GetTestDataFilePath("DiffusionImaging/Fiberfox/distortions2.dwi") ); } }; MITK_TEST_SUITE_REGISTRATION(mitkFiberfoxAddArtifactsToDwi)