diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkDftImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkDftImageFilter.cpp index 2bb093a428..8067eb1722 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkDftImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkDftImageFilter.cpp @@ -1,90 +1,111 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkDftImageFilter_txx #define __itkDftImageFilter_txx #include #include #include #include "itkDftImageFilter.h" #include #include #include #define _USE_MATH_DEFINES #include namespace itk { template< class TPixelType > DftImageFilter< TPixelType > ::DftImageFilter() { this->SetNumberOfRequiredInputs( 1 ); } +template< class TPixelType > +void DftImageFilter< TPixelType > +::BeforeThreadedGenerateData() +{ + typename InputImageType::Pointer inputImage = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); + + m_PhaseImage = OutputImageType::New(); + m_PhaseImage->SetSpacing( inputImage->GetSpacing() ); + m_PhaseImage->SetOrigin( inputImage->GetOrigin() ); + m_PhaseImage->SetDirection( inputImage->GetDirection() ); + m_PhaseImage->SetLargestPossibleRegion( inputImage->GetLargestPossibleRegion() ); + m_PhaseImage->SetBufferedRegion( inputImage->GetLargestPossibleRegion() ); + m_PhaseImage->SetRequestedRegion( inputImage->GetLargestPossibleRegion() ); + m_PhaseImage->Allocate(); + m_PhaseImage->FillBuffer(0.0); +} + template< class TPixelType > void DftImageFilter< TPixelType > ::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; typename InputImageType::Pointer inputImage = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); int szx = outputImage->GetLargestPossibleRegion().GetSize(0); int szy = outputImage->GetLargestPossibleRegion().GetSize(1); while( !oit.IsAtEnd() ) { int kx = oit.GetIndex()[0]; int ky = oit.GetIndex()[1]; if( kx < szx/2 ) kx = kx + szx/2; else kx = kx - szx/2; if( ky < szy/2 ) ky = ky + szy/2; else ky = ky - szy/2; vcl_complex s(0,0); InputIteratorType it(inputImage, inputImage->GetLargestPossibleRegion() ); while( !it.IsAtEnd() ) { int x = it.GetIndex()[0]; int y = it.GetIndex()[1]; vcl_complex f(it.Get().real(), it.Get().imag()); + s += f * exp( std::complex(0, -2 * M_PI * (kx*(double)x/szx + ky*(double)y/szy) ) ); ++it; } double magn = sqrt(s.real()*s.real()+s.imag()*s.imag()); - oit.Set(magn); + if (s.real() != 0 ) + m_PhaseImage->SetPixel(oit.GetIndex(), atan(s.imag()/s.real())); + + oit.Set(magn); ++oit; } } } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkDftImageFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkDftImageFilter.h index 6062dc8cc0..2ce218db06 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkDftImageFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkDftImageFilter.h @@ -1,75 +1,83 @@ /*=================================================================== 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 __itkDftImageFilter_h_ #define __itkDftImageFilter_h_ #include #include #include #include +#include namespace itk{ /** * \brief 2D Discrete Fourier Transform Filter (complex to real). Special issue for Fiberfox -> rearranges slice. */ template< class TPixelType > class DftImageFilter : public ImageToImageFilter< Image< vcl_complex< TPixelType > >, Image< TPixelType > > { public: typedef DftImageFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< Image< vcl_complex< TPixelType > >, Image< TPixelType > > Superclass; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(DftImageFilter, ImageToImageFilter) typedef typename Superclass::InputImageType InputImageType; typedef typename Superclass::OutputImageType OutputImageType; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; + itkGetMacro( PhaseImage, typename OutputImageType::Pointer ) + + void SetParameters( FiberfoxParameters param ){ m_Parameters = param; } + protected: DftImageFilter(); ~DftImageFilter() {} + void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType threadId); private: + typename OutputImageType::Pointer m_PhaseImage; + FiberfoxParameters m_Parameters; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkDftImageFilter.cpp" #endif #endif //__itkDftImageFilter_h_ diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.cpp index cfcbf0e488..434fecbab9 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.cpp @@ -1,225 +1,279 @@ /*=================================================================== 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_FrequencyMapSlice(NULL) , m_Z(0) , m_UseConstantRandSeed(false) , m_SpikesPerSlice(0) , m_IsBaseline(true) { 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(); + m_Spike = vcl_complex(0,0); + 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(); + outputImage->FillBuffer(m_Spike); + + m_KSpaceImage = InputImageType::New(); + m_KSpaceImage->SetLargestPossibleRegion( region ); + m_KSpaceImage->SetBufferedRegion( region ); + m_KSpaceImage->SetRequestedRegion( region ); + m_KSpaceImage->Allocate(); + m_KSpaceImage->FillBuffer(0.0); double gamma = 42576000; // Gyromagnetic ratio in Hz/T if ( m_Parameters.m_SignalGen.m_EddyStrength>0 && m_DiffusionGradientDirection.GetNorm()>0.001) { m_DiffusionGradientDirection.Normalize(); m_DiffusionGradientDirection = m_DiffusionGradientDirection * m_Parameters.m_SignalGen.m_EddyStrength/1000 * gamma; m_IsBaseline = false; } this->SetNthOutput(0, outputImage); - m_Spike = vcl_complex(0,0); m_Transform = m_Parameters.m_SignalGen.m_ImageDirection; for (int i=0; i<3; i++) for (int j=0; j<3; j++) m_Transform[i][j] *= m_Parameters.m_SignalGen.m_ImageSpacing[j]; } template< class TPixelType > void KspaceImageFilter< TPixelType > ::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_Parameters.m_SignalGen.m_tLine/kxMax; - double fromMaxEcho = - m_Parameters.m_SignalGen.m_tLine*kyMax/2; + double fromMaxEcho = -m_Parameters.m_SignalGen.m_tLine*kyMax/2-dt*kxMax/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]; + bool pf = false; + if (kIdx[1]>kyMax*m_Parameters.m_SignalGen.m_PartialFourier) + pf = true; + // dephasing time double t = 0; if (!m_Parameters.m_SignalGen.m_ReversePhase) t = fromMaxEcho + ((double)kIdx[1]*kxMax+(double)kIdx[0])*dt; else t = fromMaxEcho + ((double)(kyMax-1-kIdx[1])*kxMax+(double)kIdx[0])*dt; // 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_Parameters.m_SignalGen.m_EddyStrength>0) eddyDecay = exp(-( m_Parameters.m_SignalGen.m_tEcho/2 + t)/ m_Parameters.m_SignalGen.m_Tau) * t/1000; // calcualte signal relaxation factors std::vector< double > relaxFactor; if ( m_Parameters.m_SignalGen.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() ) + if (!pf) { - double x = it.GetIndex()[0]-xMax/2; - double y = it.GetIndex()[1]-yMax/2; + // add gibbs ringing offset (cropps k-space) + if (kx>=kxMax/2) + kx += xRingingOffset; + if (ky>=kyMax/2) + ky += yRingingOffset; - 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_Parameters.m_SignalGen.m_SignalScale, 0); - else - f += std::complex( m_CompartmentImages.at(i)->GetPixel(it.GetIndex()) * m_Parameters.m_SignalGen.m_SignalScale ); - - // simulate eddy currents and other distortions - double omega_t = 0; - if ( m_Parameters.m_SignalGen.m_EddyStrength>0 && !m_IsBaseline) + vcl_complex s(0,0); + InputIteratorType it(m_CompartmentImages.at(0), m_CompartmentImages.at(0)->GetLargestPossibleRegion() ); + while( !it.IsAtEnd() ) { - itk::Vector< double, 3 > pos; pos[0] = x; pos[1] = y; pos[2] = m_Z; - 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; + 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_Parameters.m_SignalGen.m_SignalScale, 0); + else + f += std::complex( m_CompartmentImages.at(i)->GetPixel(it.GetIndex()) * m_Parameters.m_SignalGen.m_SignalScale ); + + // simulate eddy currents and other distortions + double omega_t = 0; + if ( m_Parameters.m_SignalGen.m_EddyStrength>0 && !m_IsBaseline) + { + itk::Vector< double, 3 > pos; pos[0] = x; pos[1] = y; pos[2] = m_Z; + 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_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; } - if (m_FrequencyMapSlice.IsNotNull()) // simulate distortions - omega_t += m_FrequencyMapSlice->GetPixel(it.GetIndex())*t/1000; + s /= numPix; - if (y<-yMaxFov/2) - y += yMaxFov; - else if (y>=yMaxFov/2) - y -= yMaxFov; + 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; - // actual DFT term - s += f * exp( std::complex(0, 2 * M_PI * (kx*x/xMax + ky*y/yMaxFov + omega_t )) ); + s = vcl_complex(s.real()+m_RandGen->GetNormalVariate(0,m_Parameters.m_SignalGen.m_NoiseVariance), s.imag()+m_RandGen->GetNormalVariate(0,m_Parameters.m_SignalGen.m_NoiseVariance)); - ++it; + outputImage->SetPixel(kIdx, s); + m_KSpaceImage->SetPixel(oit.GetIndex(), sqrt(s.imag()*s.imag()+s.real()*s.real()) ); } - s /= numPix; - 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 + ImageRegionIterator< OutputImageType > oit(outputImage, outputImage->GetLargestPossibleRegion()); + while( !oit.IsAtEnd() ) // use hermitian k-space symmetry to fill empty k-space parts resulting from partial fourier acquisition + { + itk::Index< 2 > kIdx; + kIdx[0] = oit.GetIndex()[0]; + kIdx[1] = oit.GetIndex()[1]; + + if (kIdx[1]>kyMax*m_Parameters.m_SignalGen.m_PartialFourier) + { + 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; + + if (oit.GetIndex()[1]%2 == 1) // reverse readout direction and add ghosting + kIdx[0] = kxMax-kIdx[0]-1; // reverse readout direction + + itk::Index< 2 > kIdx2; + kIdx2[0] = (int)(kxMax-kIdx[0])%(int)kxMax; + kIdx2[1] = (int)(kyMax-kIdx[1])%(int)kyMax; + + vcl_complex s = outputImage->GetPixel(kIdx2); + s = vcl_complex(s.real(), -s.imag()); + outputImage->SetPixel(kIdx, s); + m_KSpaceImage->SetPixel(oit.GetIndex(), sqrt(s.imag()*s.imag()+s.real()*s.real()) ); + } + ++oit; + } + m_Spike *= m_Parameters.m_SignalGen.m_SpikeAmplitude; itk::Index< 2 > spikeIdx; for (unsigned 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 1e2ad4c9ff..6ef1bc29d2 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.h @@ -1,124 +1,127 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ /*=================================================================== This file is based heavily on a corresponding ITK filter. ===================================================================*/ #ifndef __itkKspaceImageFilter_h_ #define __itkKspaceImageFilter_h_ #include #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. */ itkFactorylessNewMacro(Self) itkCloneMacro(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( 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( OutSize, itk::Size<2> ) ///< Output slice size. Can be different from input size, e.g. if Gibbs ringing is enabled. itkSetMacro( UseConstantRandSeed, bool ) ///< Use constant seed for random generator for reproducible results. + itkGetMacro( KSpaceImage, typename InputImageType::Pointer ) ///< k-space magnitude image 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 SetT1( std::vector< double > t1Vector ) { m_T1=t1Vector; } ///< One T1 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); void AfterThreadedGenerateData(); FiberfoxParameters m_Parameters; typename InputImageType::Pointer m_FrequencyMapSlice; vector< double > m_T2; vector< double > m_T1; vector< InputImagePointerType > m_CompartmentImages; itk::Vector m_DiffusionGradientDirection; double m_Z; 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; + typename InputImageType::Pointer m_KSpaceImage; + 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 7906df9dcb..a7dfe30e9a 100755 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp @@ -1,1325 +1,1359 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "itkTractsToDWIImageFilter.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include namespace itk { template< class PixelType > TractsToDWIImageFilter< PixelType >::TractsToDWIImageFilter() : m_FiberBundle(NULL) , m_StatusText("") , m_UseConstantRandSeed(false) , m_RandGen(itk::Statistics::MersenneTwisterRandomVariateGenerator::New()) { m_RandGen->SetSeed(); } template< class PixelType > TractsToDWIImageFilter< PixelType >::~TractsToDWIImageFilter() { } template< class PixelType > TractsToDWIImageFilter< PixelType >::DoubleDwiType::Pointer TractsToDWIImageFilter< PixelType >::DoKspaceStuff( std::vector< DoubleDwiType::Pointer >& images ) { int numFiberCompartments = m_Parameters.m_FiberModelList.size(); // create slice object ImageRegion<2> sliceRegion; sliceRegion.SetSize(0, m_WorkingImageRegion.GetSize()[0]); sliceRegion.SetSize(1, m_WorkingImageRegion.GetSize()[1]); Vector< double, 2 > sliceSpacing; sliceSpacing[0] = m_WorkingSpacing[0]; sliceSpacing[1] = m_WorkingSpacing[1]; // frequency map slice SliceType::Pointer fMapSlice = NULL; if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull()) { fMapSlice = SliceType::New(); ImageRegion<2> region; region.SetSize(0, m_WorkingImageRegion.GetSize()[0]); region.SetSize(1, m_WorkingImageRegion.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_SignalGen.m_ImageSpacing ); newImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin ); newImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); newImage->SetLargestPossibleRegion( m_CroppedRegion ); newImage->SetBufferedRegion( m_CroppedRegion ); newImage->SetRequestedRegion( m_CroppedRegion ); newImage->SetVectorLength( images.at(0)->GetVectorLength() ); newImage->Allocate(); + 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_CroppedRegion ); + m_PhaseImage->SetBufferedRegion( m_CroppedRegion ); + m_PhaseImage->SetRequestedRegion( m_CroppedRegion ); + m_PhaseImage->SetVectorLength( images.at(0)->GetVectorLength() ); + m_PhaseImage->Allocate(); + + 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_CroppedRegion ); + m_KspaceImage->SetBufferedRegion( m_CroppedRegion ); + m_KspaceImage->SetRequestedRegion( m_CroppedRegion ); + m_KspaceImage->SetVectorLength( images.at(0)->GetVectorLength() ); + m_KspaceImage->Allocate(); + std::vector< unsigned int > spikeVolume; for (unsigned int i=0; iGetIntegerVariate()%(images.at(0)->GetVectorLength())); std::sort (spikeVolume.begin(), spikeVolume.end()); std::reverse (spikeVolume.begin(), spikeVolume.end()); m_StatusText += "0% 10 20 30 40 50 60 70 80 90 100%\n"; m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*"; unsigned long lastTick = 0; boost::progress_display disp(2*images.at(0)->GetVectorLength()*images.at(0)->GetLargestPossibleRegion().GetSize(2)); for (unsigned int g=0; gGetVectorLength(); g++) { std::vector< unsigned int > spikeSlice; while (!spikeVolume.empty() && spikeVolume.back()==g) { spikeSlice.push_back(m_RandGen->GetIntegerVariate()%images.at(0)->GetLargestPossibleRegion().GetSize(2)); spikeVolume.pop_back(); } std::sort (spikeSlice.begin(), spikeSlice.end()); std::reverse (spikeSlice.begin(), spikeSlice.end()); for (unsigned int z=0; zGetLargestPossibleRegion().GetSize(2); z++) { std::vector< SliceType::Pointer > compartmentSlices; std::vector< double > t2Vector; std::vector< double > t1Vector; 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_SignalGen.m_FrequencyMap->GetPixel(index3D)); } compartmentSlices.push_back(slice); t2Vector.push_back(signalModel->GetT2()); t1Vector.push_back(signalModel->GetT1()); } if (this->GetAbortGenerateData()) return NULL; // create k-sapce (inverse fourier transform slices) itk::Size<2> outSize; outSize.SetElement(0, m_CroppedRegion.GetSize(0)); outSize.SetElement(1, m_CroppedRegion.GetSize(1)); itk::KspaceImageFilter< SliceType::PixelType >::Pointer idft = itk::KspaceImageFilter< SliceType::PixelType >::New(); idft->SetCompartmentImages(compartmentSlices); idft->SetT2(t2Vector); idft->SetT1(t1Vector); idft->SetUseConstantRandSeed(m_UseConstantRandSeed); idft->SetParameters(m_Parameters); idft->SetZ((double)z-(double)images.at(0)->GetLargestPossibleRegion().GetSize(2)/2.0); idft->SetDiffusionGradientDirection(m_Parameters.m_SignalGen.GetGradientDirection(g)); idft->SetFrequencyMapSlice(fMapSlice); idft->SetOutSize(outSize); int numSpikes = 0; while (!spikeSlice.empty() && spikeSlice.back()==z) { numSpikes++; spikeSlice.pop_back(); } idft->SetSpikesPerSlice(numSpikes); idft->Update(); ComplexSliceType::Pointer fSlice; fSlice = idft->GetOutput(); ++disp; unsigned long newTick = 50*disp.count()/disp.expected_count(); for (unsigned long tick = 0; tick<(newTick-lastTick); tick++) m_StatusText += "*"; lastTick = newTick; // fourier transform slice SliceType::Pointer newSlice; itk::DftImageFilter< SliceType::PixelType >::Pointer dft = itk::DftImageFilter< SliceType::PixelType >::New(); dft->SetInput(fSlice); dft->Update(); newSlice = dft->GetOutput(); // put slice back into channel g for (unsigned int y=0; yGetLargestPossibleRegion().GetSize(1); y++) for (unsigned int x=0; xGetLargestPossibleRegion().GetSize(0); x++) { DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z; SliceType::IndexType index2D; index2D[0]=x; index2D[1]=y; DoubleDwiType::PixelType pix3D = newImage->GetPixel(index3D); pix3D[g] = newSlice->GetPixel(index2D); newImage->SetPixel(index3D, pix3D); + + // phase image + { + DoubleDwiType::PixelType pix3D = m_PhaseImage->GetPixel(index3D); + pix3D[g] = dft->GetPhaseImage()->GetPixel(index2D); + m_PhaseImage->SetPixel(index3D, pix3D); + } + + // phase image + { + DoubleDwiType::PixelType pix3D = m_KspaceImage->GetPixel(index3D); + pix3D[g] = idft->GetKSpaceImage()->GetPixel(index2D); + m_KspaceImage->SetPixel(index3D, pix3D); + } } ++disp; newTick = 50*disp.count()/disp.expected_count(); for (unsigned long tick = 0; tick<(newTick-lastTick); tick++) m_StatusText += "*"; lastTick = newTick; } } m_StatusText += "\n\n"; return newImage; } template< class PixelType > void TractsToDWIImageFilter< PixelType >::CheckVolumeFractionImages() { // check for fiber volume fraction maps for (int i=0; iGetVolumeFractionImage().IsNotNull()) m_StatusText += "Using volume fraction map for fiber compartment " + boost::lexical_cast(i+1) + "\n"; // check for non-fiber volume fraction maps for (int i=0; iGetVolumeFractionImage().IsNotNull()) m_StatusText += "Using volume fraction map for non-fiber compartment " + boost::lexical_cast(i+1) + "\n"; // check volume fraction images if (m_Parameters.m_FiberModelList.size()>2) for (int i=0; iGetVolumeFractionImage().IsNull()) itkExceptionMacro("More than two fiber compartment selected but no volume fraction maps set!"); if (m_Parameters.m_NonFiberModelList.size()>1) for (int i=0; iGetVolumeFractionImage().IsNull()) itkExceptionMacro("More than one non-fiber compartment selected but no volume fraction maps set!"); // initialize the images that store the output volume fraction of each compartment m_VolumeFractions.clear(); for (int i=0; iSetSpacing( 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() { // initialize output dwi image m_CroppedRegion = m_Parameters.m_SignalGen.m_ImageRegion; m_CroppedRegion.SetSize(1, m_CroppedRegion.GetSize(1)*m_Parameters.m_SignalGen.m_CroppingFactor); itk::Point shiftedOrigin = m_Parameters.m_SignalGen.m_ImageOrigin; shiftedOrigin[1] += (m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1)-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_CroppedRegion ); m_OutputImage->SetBufferedRegion( m_CroppedRegion ); m_OutputImage->SetRequestedRegion( 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); if ( m_CroppedRegion.GetSize(0)%2 == 1 ) m_CroppedRegion.SetSize(0, m_CroppedRegion.GetSize(0)+1); if ( m_CroppedRegion.GetSize(1)%2 == 1 ) m_CroppedRegion.SetSize(1, m_CroppedRegion.GetSize(1)+1); // ADJUST GEOMETRY FOR FURTHER PROCESSING // is input slize size a power of two? unsigned int x=m_Parameters.m_SignalGen.m_ImageRegion.GetSize(0); unsigned int y=m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1); ItkDoubleImgType::SizeType pad; pad[0]=x%2; pad[1]=y%2; pad[2]=0; m_Parameters.m_SignalGen.m_ImageRegion.SetSize(0, x+pad[0]); m_Parameters.m_SignalGen.m_ImageRegion.SetSize(1, y+pad[1]); if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull() && (pad[0]>0 || pad[1]>0)) { itk::ConstantPadImageFilter::Pointer zeroPadder = itk::ConstantPadImageFilter::New(); zeroPadder->SetInput(m_Parameters.m_SignalGen.m_FrequencyMap); zeroPadder->SetConstant(0); zeroPadder->SetPadUpperBound(pad); zeroPadder->Update(); m_Parameters.m_SignalGen.m_FrequencyMap = zeroPadder->GetOutput(); } if (m_Parameters.m_SignalGen.m_MaskImage.IsNotNull() && (pad[0]>0 || pad[1]>0)) { itk::ConstantPadImageFilter::Pointer zeroPadder = itk::ConstantPadImageFilter::New(); zeroPadder->SetInput(m_Parameters.m_SignalGen.m_MaskImage); zeroPadder->SetConstant(0); zeroPadder->SetPadUpperBound(pad); zeroPadder->Update(); m_Parameters.m_SignalGen.m_MaskImage = zeroPadder->GetOutput(); } // 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; iSetSpacing( 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); } // resample mask image and frequency map to fit upsampled geometry if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging) { if (m_Parameters.m_SignalGen.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_SignalGen.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_SignalGen.m_MaskImage); resampler->SetSize(m_WorkingImageRegion.GetSize()); resampler->SetOutputSpacing(m_WorkingSpacing); resampler->SetOutputOrigin(m_WorkingOrigin); itk::NearestNeighborInterpolateImageFunction::Pointer nn_interpolator = itk::NearestNeighborInterpolateImageFunction::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()) { itk::ResampleImageFilter::Pointer resampler = itk::ResampleImageFilter::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); itk::NearestNeighborInterpolateImageFunction::Pointer nn_interpolator = itk::NearestNeighborInterpolateImageFunction::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 m_StatusText += "No tissue mask set\n"; MITK_INFO << "No tissue mask set"; 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 { if (m_Parameters.m_SignalGen.m_MaskImage->GetLargestPossibleRegion()!=m_Parameters.m_SignalGen.m_ImageRegion) itkExceptionMacro("Mask image and specified DWI geometry are not matching!"); m_StatusText += "Using tissue mask\n"; MITK_INFO << "Using tissue mask"; } if (m_Parameters.m_SignalGen.m_DoAddMotion) { std::string fileName = "fiberfox_motion_0.log"; std::string filePath = mitk::IOUtil::GetTempPath(); if (m_Parameters.m_Misc.m_OutputPath.size()>0) filePath = m_Parameters.m_Misc.m_OutputPath; int c = 1; while (itksys::SystemTools::FileExists((filePath+fileName).c_str())) { fileName = "fiberfox_motion_"; fileName += boost::lexical_cast(c); fileName += ".log"; c++; } m_Logfile.open((filePath+fileName).c_str()); m_Logfile << "0 rotation: 0,0,0; translation: 0,0,0\n"; if (m_Parameters.m_SignalGen.m_DoRandomizeMotion) { m_StatusText += "Adding random motion artifacts:\n"; m_StatusText += "Maximum rotation: +/-" + boost::lexical_cast(m_Parameters.m_SignalGen.m_Rotation) + "°\n"; m_StatusText += "Maximum translation: +/-" + boost::lexical_cast(m_Parameters.m_SignalGen.m_Translation) + "mm\n"; } else { m_StatusText += "Adding linear motion artifacts:\n"; m_StatusText += "Maximum rotation: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_Rotation) + "°\n"; m_StatusText += "Maximum translation: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_Translation) + "mm\n"; } m_StatusText += "Motion logfile: " + (filePath+fileName) + "\n"; MITK_INFO << "Adding motion artifacts"; MITK_INFO << "Maximum rotation: " << m_Parameters.m_SignalGen.m_Rotation; MITK_INFO << "Maxmimum translation: " << m_Parameters.m_SignalGen.m_Translation; } // get transform for motion artifacts m_Rotation = m_Parameters.m_SignalGen.m_Rotation/m_Parameters.m_SignalGen.GetNumVolumes(); m_Translation = m_Parameters.m_SignalGen.m_Translation/m_Parameters.m_SignalGen.GetNumVolumes(); // creat image to hold transformed mask (motion artifact) m_TransformedMaskImage = ItkUcharImgType::New(); itk::ImageDuplicator::Pointer duplicator = itk::ImageDuplicator::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 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(); itk::ResampleImageFilter::Pointer upsampler = itk::ResampleImageFilter::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); itk::NearestNeighborInterpolateImageFunction::Pointer nn_interpolator = itk::NearestNeighborInterpolateImageFunction::New(); upsampler->SetInterpolator(nn_interpolator); upsampler->Update(); m_UpsampledMaskImage = upsampler->GetOutput(); } template< class PixelType > void TractsToDWIImageFilter< PixelType >::InitializeFiberData() { // resample fiber bundle for sufficient voxel coverage m_StatusText += "\n"+this->GetTime()+" > Resampling fibers ...\n"; segmentVolume = 0.0001; float minSpacing = 1; if(m_WorkingSpacing[0]GetDeepCopy(); // working copy is needed because we need to resample the fibers but do not want to change the original bundle double volumeAccuracy = 10; m_FiberBundleWorkingCopy->ResampleSpline(minSpacing/volumeAccuracy); mmRadius = m_Parameters.m_SignalGen.m_AxonRadius/1000; if (mmRadius>0) segmentVolume = M_PI*mmRadius*mmRadius*minSpacing/volumeAccuracy; m_FiberBundleTransformed = m_FiberBundleWorkingCopy; // a secon fiber bundle is needed to store the transformed version of the m_FiberBundleWorkingCopy } template< class PixelType > void TractsToDWIImageFilter< PixelType >::GenerateData() { m_TimeProbe.Start(); m_StatusText = "Starting simulation\n"; // check input data if (m_FiberBundle.IsNull()) itkExceptionMacro("Input fiber bundle is NULL!"); if (m_Parameters.m_FiberModelList.empty()) itkExceptionMacro("No diffusion model for fiber compartments defined! At least one fiber compartment is necessary."); if (m_Parameters.m_NonFiberModelList.empty()) itkExceptionMacro("No diffusion model for non-fiber compartments defined! At least one non-fiber compartment is necessary."); 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(); 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(); m_StatusText += "\n"+this->GetTime()+" > Generating " + boost::lexical_cast(numFiberCompartments+numNonFiberCompartments) + "-compartment diffusion-weighted signal.\n"; MITK_INFO << "Generating " << numFiberCompartments+numNonFiberCompartments << "-compartment diffusion-weighted signal."; int numFibers = m_FiberBundleWorkingCopy->GetNumFibers(); boost::progress_display disp(numFibers*m_Parameters.m_SignalGen.GetNumVolumes()); switch (m_Parameters.m_SignalGen.m_DiffusionDirectionMode) { case(SignalGenerationParameters::FIBER_TANGENT_DIRECTIONS): // use fiber tangent directions to determine diffusion direction { m_StatusText += "0% 10 20 30 40 50 60 70 80 90 100%\n"; m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*"; for (unsigned int g=0; gSetSeed(signalModelSeed); for (int i=0; iSetSeed(signalModelSeed); // storing voxel-wise intra-axonal volume in mm³ ItkDoubleImgType::Pointer 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; vtkPolyData* fiberPolyData = m_FiberBundleTransformed->GetFiberPolyData(); // generate fiber signal (if there are any fiber models present) if (!m_Parameters.m_FiberModelList.empty()) for( int i=0; iGetFiberWeight(i); vtkCell* cell = fiberPolyData->GetCell(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; m_TransformedMaskImage->TransformPhysicalPointToIndex(vertex, idx); m_TransformedMaskImage->TransformPhysicalPointToContinuousIndex(vertex, contIndex); if (!m_TransformedMaskImage->GetLargestPossibleRegion().IsInside(idx) || m_TransformedMaskImage->GetPixel(idx)<=0) continue; // generate signal for each fiber compartment for (int k=0; kSetFiberDirection(dir); DoubleDwiType::PixelType pix = m_CompartmentImages.at(k)->GetPixel(idx); pix[g] += fiberWeight*segmentVolume*m_Parameters.m_FiberModelList[k]->SimulateMeasurement(g); m_CompartmentImages.at(k)->SetPixel(idx, pix); } // update fiber volume image double vol = intraAxonalVolumeImage->GetPixel(idx) + segmentVolume*fiberWeight; intraAxonalVolumeImage->SetPixel(idx, vol); if (vol>maxVolume) // we assume that the first volume is always unweighted! maxVolume = vol; } // progress report ++disp; unsigned long newTick = 50*disp.count()/disp.expected_count(); for (unsigned int tick = 0; tick<(newTick-lastTick); tick++) m_StatusText += "*"; lastTick = newTick; } // generate non-fiber signal ImageRegionIterator it3(m_TransformedMaskImage, m_TransformedMaskImage->GetLargestPossibleRegion()); double fact = 1; // density correction factor in mm³ if (m_Parameters.m_SignalGen.m_AxonRadius<0.0001 || maxVolume>m_VoxelVolume) // the fullest voxel is always completely full fact = m_VoxelVolume/maxVolume; while(!it3.IsAtEnd()) { if (it3.Get()>0) { DoubleDwiType::IndexType index = it3.GetIndex(); itk::Point point; m_TransformedMaskImage->TransformIndexToPhysicalPoint(index, point); if (m_Parameters.m_SignalGen.m_DoAddMotion) { if (m_Parameters.m_SignalGen.m_DoRandomizeMotion && g>0) point = m_FiberBundleWorkingCopy->TransformPoint(point.GetVnlVector(), -m_Rotation[0],-m_Rotation[1],-m_Rotation[2],-m_Translation[0],-m_Translation[1],-m_Translation[2]); else if (g>=0) point = m_FiberBundleWorkingCopy->TransformPoint(point.GetVnlVector(), -m_Rotation[0]*g,-m_Rotation[1]*g,-m_Rotation[2]*g,-m_Translation[0]*g,-m_Translation[1]*g,-m_Translation[2]*g); } double iAxVolume = intraAxonalVolumeImage->GetPixel(index); // if volume fraction image is set use it, otherwise use scaling factor to obtain one full fiber voxel double fact2 = fact; if (m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()!=nullptr && iAxVolume>0.0001) { DoubleDwiType::IndexType newIndex; m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()->TransformPhysicalPointToIndex(point, newIndex); if (m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()->GetLargestPossibleRegion().IsInside(newIndex)) fact2 = m_VoxelVolume*m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()->GetPixel(newIndex)/iAxVolume; } // adjust intra-axonal image value for (int i=0; iGetPixel(index); pix[g] *= fact2; m_CompartmentImages.at(i)->SetPixel(index, pix); } // simulate other compartments SimulateExtraAxonalSignal(index, iAxVolume*fact2, g); } ++it3; } // move fibers SimulateMotion(g); } break; } case (SignalGenerationParameters::MAIN_FIBER_DIRECTIONS): // use main fiber directions to determine voxel-wise diffusion directions { typedef itk::Image< itk::Vector< float, 3>, 3 > ItkDirectionImage3DType; typedef itk::VectorContainer< unsigned int, ItkDirectionImage3DType::Pointer > ItkDirectionImageContainerType; // calculate main fiber directions itk::TractsToVectorImageFilter::Pointer fOdfFilter = itk::TractsToVectorImageFilter::New(); fOdfFilter->SetFiberBundle(m_FiberBundleTransformed); fOdfFilter->SetMaskImage(m_TransformedMaskImage); fOdfFilter->SetAngularThreshold(cos(m_Parameters.m_SignalGen.m_FiberSeparationThreshold*M_PI/180.0)); fOdfFilter->SetNormalizeVectors(false); fOdfFilter->SetUseWorkingCopy(true); fOdfFilter->SetSizeThreshold(0); fOdfFilter->SetMaxNumDirections(3); fOdfFilter->Update(); ItkDirectionImageContainerType::Pointer directionImageContainer = fOdfFilter->GetDirectionImageContainer(); // allocate image storing intra-axonal volume fraction information ItkDoubleImgType::Pointer 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); // determine intra-axonal volume fraction using the tract density itk::TractDensityImageFilter< ItkDoubleImgType >::Pointer tdiFilter = itk::TractDensityImageFilter< ItkDoubleImgType >::New(); tdiFilter->SetFiberBundle(m_FiberBundleTransformed); tdiFilter->SetBinaryOutput(false); tdiFilter->SetOutputAbsoluteValues(false); tdiFilter->SetInputImage(intraAxonalVolumeImage); tdiFilter->SetUseImageGeometry(true); tdiFilter->Update(); intraAxonalVolumeImage = tdiFilter->GetOutput(); m_StatusText += "0% 10 20 30 40 50 60 70 80 90 100%\n"; m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*"; boost::progress_display disp(m_TransformedMaskImage->GetLargestPossibleRegion().GetNumberOfPixels()*m_Parameters.m_SignalGen.GetNumVolumes()); for (unsigned int g=0; gSetSeed(signalModelSeed); for (int i=0; iSetSeed(signalModelSeed); if (m_Parameters.m_SignalGen.m_DoAddMotion && g>0) // if fibers have moved we need a new TDI and new directions { fOdfFilter->SetFiberBundle(m_FiberBundleTransformed); fOdfFilter->SetMaskImage(m_TransformedMaskImage); fOdfFilter->Update(); directionImageContainer = fOdfFilter->GetDirectionImageContainer(); tdiFilter->SetFiberBundle(m_FiberBundleTransformed); tdiFilter->Update(); intraAxonalVolumeImage = tdiFilter->GetOutput(); } ImageRegionIterator< ItkUcharImgType > it(m_TransformedMaskImage, m_TransformedMaskImage->GetLargestPossibleRegion()); while(!it.IsAtEnd()) { ++disp; unsigned long newTick = 50*disp.count()/disp.expected_count(); for (unsigned int tick = 0; tick<(newTick-lastTick); tick++) m_StatusText += "*"; lastTick = newTick; if (this->GetAbortGenerateData()) { m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n"; return; } if (it.Get()>0) { // generate fiber signal for (int c=0; cGetPixel(it.GetIndex()); for (unsigned int i=0; iSize(); i++) { itk::Vector< double, 3> dir; dir.CastFrom(directionImageContainer->GetElement(i)->GetPixel(it.GetIndex())); double norm = dir.GetNorm(); if (norm>0.0001) { m_Parameters.m_FiberModelList.at(c)->SetFiberDirection(dir); pix[g] += m_Parameters.m_FiberModelList.at(c)->SimulateMeasurement(g)*norm; count++; } } if (count>0) pix[g] /= count; pix[g] *= intraAxonalVolumeImage->GetPixel(it.GetIndex())*m_VoxelVolume; m_CompartmentImages.at(c)->SetPixel(it.GetIndex(), pix); } // simulate other compartments SimulateExtraAxonalSignal(it.GetIndex(), intraAxonalVolumeImage->GetPixel(it.GetIndex())*m_VoxelVolume, g); } ++it; } SimulateMotion(g); } itk::ImageFileWriter< ItkUcharImgType >::Pointer wr = itk::ImageFileWriter< ItkUcharImgType >::New(); wr->SetInput(fOdfFilter->GetNumDirectionsImage()); wr->SetFileName(mitk::IOUtil::GetTempPath()+"/NumDirections_MainFiberDirections.nrrd"); wr->Update(); break; } case (SignalGenerationParameters::RANDOM_DIRECTIONS): { ItkUcharImgType::Pointer numDirectionsImage = ItkUcharImgType::New(); numDirectionsImage->SetSpacing( m_WorkingSpacing ); numDirectionsImage->SetOrigin( m_WorkingOrigin ); numDirectionsImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); numDirectionsImage->SetLargestPossibleRegion( m_WorkingImageRegion ); numDirectionsImage->SetBufferedRegion( m_WorkingImageRegion ); numDirectionsImage->SetRequestedRegion( m_WorkingImageRegion ); numDirectionsImage->Allocate(); numDirectionsImage->FillBuffer(0); double sepAngle = cos(m_Parameters.m_SignalGen.m_FiberSeparationThreshold*M_PI/180.0); m_StatusText += "0% 10 20 30 40 50 60 70 80 90 100%\n"; m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*"; boost::progress_display disp(m_TransformedMaskImage->GetLargestPossibleRegion().GetNumberOfPixels()); ImageRegionIterator it(m_TransformedMaskImage, m_TransformedMaskImage->GetLargestPossibleRegion()); while(!it.IsAtEnd()) { ++disp; unsigned long newTick = 50*disp.count()/disp.expected_count(); for (unsigned int tick = 0; tick<(newTick-lastTick); tick++) m_StatusText += "*"; lastTick = newTick; if (this->GetAbortGenerateData()) { m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n"; return; } if (it.Get()>0) { int numFibs = m_RandGen->GetIntegerVariate(2)+1; DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(it.GetIndex()); double volume = m_RandGen->GetVariateWithClosedRange(0.3); // double sum = 0; std::vector< double > fractions; for (int i=0; iGetVariateWithClosedRange(0.5)); // sum += fractions.at(i); } // for (int i=0; i > directions; for (int i=0; iGetVariateWithClosedRange(2)-1.0; fib[1] = m_RandGen->GetVariateWithClosedRange(2)-1.0; fib[2] = m_RandGen->GetVariateWithClosedRange(2)-1.0; fib.Normalize(); double min = 0; for (unsigned int d=0; dmin) min = angle; } if (minSetFiberDirection(fib); pix += m_Parameters.m_FiberModelList.at(0)->SimulateMeasurement()*fractions[i]; directions.push_back(fib); } else i--; } pix *= (1-volume); m_CompartmentImages.at(0)->SetPixel(it.GetIndex(), pix); // CSF/GM { pix += volume*m_Parameters.m_NonFiberModelList.at(0)->SimulateMeasurement(); } numDirectionsImage->SetPixel(it.GetIndex(), numFibs); } ++it; } itk::ImageFileWriter< ItkUcharImgType >::Pointer wr = itk::ImageFileWriter< ItkUcharImgType >::New(); wr->SetInput(numDirectionsImage); wr->SetFileName(mitk::IOUtil::GetTempPath()+"/NumDirections_RandomDirections.nrrd"); wr->Update(); } } if (m_Logfile.is_open()) { m_Logfile << "DONE"; m_Logfile.close(); } m_StatusText += "\n\n"; if (this->GetAbortGenerateData()) { m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n"; return; } DoubleDwiType::Pointer doubleOutImage; if ( m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition ) // do k-space stuff { m_StatusText += this->GetTime()+" > Adjusting complex signal\n"; MITK_INFO << "Adjusting complex signal:"; if (m_Parameters.m_SignalGen.m_DoSimulateRelaxation) m_StatusText += "Simulating signal relaxation\n"; if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull()) m_StatusText += "Simulating distortions\n"; if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging) m_StatusText += "Simulating ringing artifacts\n"; if (m_Parameters.m_SignalGen.m_EddyStrength>0) m_StatusText += "Simulating eddy currents\n"; if (m_Parameters.m_SignalGen.m_Spikes>0) m_StatusText += "Simulating spikes\n"; if (m_Parameters.m_SignalGen.m_CroppingFactor<1.0) m_StatusText += "Simulating aliasing artifacts\n"; if (m_Parameters.m_SignalGen.m_KspaceLineOffset>0) m_StatusText += "Simulating ghosts\n"; doubleOutImage = DoKspaceStuff(m_CompartmentImages); m_Parameters.m_SignalGen.m_SignalScale = 1; // already scaled in DoKspaceStuff } else // don't do k-space stuff, just sum compartments { m_StatusText += this->GetTime()+" > Summing compartments\n"; MITK_INFO << "Summing compartments"; doubleOutImage = m_CompartmentImages.at(0); for (unsigned int i=1; i::Pointer adder = itk::AddImageFilter< DoubleDwiType, DoubleDwiType, DoubleDwiType>::New(); adder->SetInput1(doubleOutImage); adder->SetInput2(m_CompartmentImages.at(i)); adder->Update(); doubleOutImage = adder->GetOutput(); } } if (this->GetAbortGenerateData()) { m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n"; return; } m_StatusText += this->GetTime()+" > Finalizing image\n"; MITK_INFO << "Finalizing image"; if (m_Parameters.m_SignalGen.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 (m_OutputImage, m_OutputImage->GetLargestPossibleRegion()); DoubleDwiType::PixelType signal; signal.SetSize(m_Parameters.m_SignalGen.GetNumVolumes()); boost::progress_display disp2(m_OutputImage->GetLargestPossibleRegion().GetNumberOfPixels()); m_StatusText += "0% 10 20 30 40 50 60 70 80 90 100%\n"; m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*"; lastTick = 0; while(!it4.IsAtEnd()) { if (this->GetAbortGenerateData()) { m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n"; return; } ++disp2; unsigned long newTick = 50*disp2.count()/disp2.expected_count(); for (unsigned long tick = 0; tick<(newTick-lastTick); tick++) m_StatusText += "*"; lastTick = newTick; typename OutputImageType::IndexType index = it4.GetIndex(); signal = doubleOutImage->GetPixel(index)*m_Parameters.m_SignalGen.m_SignalScale; if (m_Parameters.m_NoiseModel!=NULL) m_Parameters.m_NoiseModel->AddNoise(signal); for (unsigned int i=0; i0) signal[i] = floor(signal[i]+0.5); else signal[i] = ceil(signal[i]-0.5); if ( (!m_Parameters.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]SetNthOutput(0, m_OutputImage); m_StatusText += "\n\n"; m_StatusText += "Finished simulation\n"; m_StatusText += "Simulation time: "+GetTime(); m_TimeProbe.Stop(); } template< class PixelType > void TractsToDWIImageFilter< PixelType >::SimulateMotion(int g) { if (m_Parameters.m_SignalGen.m_DoAddMotion && gGetDeepCopy(); 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]; } // rotate mask image if (m_MaskImageSet) { ImageRegionIterator maskIt(m_UpsampledMaskImage, m_UpsampledMaskImage->GetLargestPossibleRegion()); m_TransformedMaskImage->FillBuffer(0); while(!maskIt.IsAtEnd()) { if (maskIt.Get()<=0) { ++maskIt; continue; } DoubleDwiType::IndexType index = maskIt.GetIndex(); itk::Point point; m_UpsampledMaskImage->TransformIndexToPhysicalPoint(index, point); if (m_Parameters.m_SignalGen.m_DoRandomizeMotion) point = m_FiberBundleWorkingCopy->TransformPoint(point.GetVnlVector(), m_Rotation[0],m_Rotation[1],m_Rotation[2],m_Translation[0],m_Translation[1],m_Translation[2]); else point = m_FiberBundleWorkingCopy->TransformPoint(point.GetVnlVector(), m_Rotation[0]*(g+1),m_Rotation[1]*(g+1),m_Rotation[2]*(g+1),m_Translation[0]*(g+1),m_Translation[1]*(g+1),m_Translation[2]*(g+1)); m_TransformedMaskImage->TransformPhysicalPointToIndex(point, index); if (m_TransformedMaskImage->GetLargestPossibleRegion().IsInside(index)) m_TransformedMaskImage->SetPixel(index,100); ++maskIt; } } // rotate fibers if (m_Logfile.is_open()) { m_Logfile << g+1 << " rotation: " << m_Rotation[0] << "," << m_Rotation[1] << "," << m_Rotation[2] << ";"; m_Logfile << " translation: " << m_Translation[0] << "," << m_Translation[1] << "," << m_Translation[2] << "\n"; } m_FiberBundleTransformed->TransformFibers(m_Rotation[0],m_Rotation[1],m_Rotation[2],m_Translation[0],m_Translation[1],m_Translation[2]); } } template< class PixelType > void TractsToDWIImageFilter< PixelType >::SimulateExtraAxonalSignal(ItkUcharImgType::IndexType index, double intraAxonalVolume, int g) { int numFiberCompartments = m_Parameters.m_FiberModelList.size(); int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size(); if (intraAxonalVolume>0.0001 && m_Parameters.m_SignalGen.m_DoDisablePartialVolume) // only fiber in voxel { DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(index); if (g>=0) pix[g] *= m_VoxelVolume/intraAxonalVolume; else pix *= m_VoxelVolume/intraAxonalVolume; m_CompartmentImages.at(0)->SetPixel(index, pix); if (g==0) m_VolumeFractions.at(0)->SetPixel(index, 1); for (int i=1; iGetPixel(index); if (g>=0) pix[g] = 0.0; else pix.Fill(0.0); m_CompartmentImages.at(i)->SetPixel(index, pix); } } else { if (g==0) m_VolumeFractions.at(0)->SetPixel(index, intraAxonalVolume/m_VoxelVolume); // get non-transformed point (remove headmotion tranformation) // this point can then be transformed to each of the original images, regardless of their geometry itk::Point point; m_TransformedMaskImage->TransformIndexToPhysicalPoint(index, point); if (m_Parameters.m_SignalGen.m_DoAddMotion) { if (m_Parameters.m_SignalGen.m_DoRandomizeMotion && g>0) point = m_FiberBundleWorkingCopy->TransformPoint(point.GetVnlVector(), -m_Rotation[0],-m_Rotation[1],-m_Rotation[2],-m_Translation[0],-m_Translation[1],-m_Translation[2]); else if (g>=0) point = m_FiberBundleWorkingCopy->TransformPoint(point.GetVnlVector(), -m_Rotation[0]*g,-m_Rotation[1]*g,-m_Rotation[2]*g,-m_Translation[0]*g,-m_Translation[1]*g,-m_Translation[2]*g); } if (m_Parameters.m_SignalGen.m_DoDisablePartialVolume) { int maxVolumeIndex = 0; double maxWeight = 0; for (int i=0; i1) { // ItkUcharImgType::IndexType maskIndex; // m_UpsampledMaskImage->TransformPhysicalPointToIndex(point, maskIndex); // if (!m_UpsampledMaskImage->GetLargestPossibleRegion().IsInside(maskIndex) || m_UpsampledMaskImage->GetPixel(maskIndex)==0) // continue; DoubleDwiType::IndexType newIndex; m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->TransformPhysicalPointToIndex(point, newIndex); if (!m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetLargestPossibleRegion().IsInside(newIndex)) continue; weight = m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetPixel(newIndex); } if (weight>maxWeight) { maxWeight = weight; maxVolumeIndex = i; } } DoubleDwiType::Pointer doubleDwi = m_CompartmentImages.at(maxVolumeIndex+numFiberCompartments); DoubleDwiType::PixelType pix = doubleDwi->GetPixel(index); if (g>=0) pix[g] += m_Parameters.m_NonFiberModelList[maxVolumeIndex]->SimulateMeasurement(g); else pix += m_Parameters.m_NonFiberModelList[maxVolumeIndex]->SimulateMeasurement(); doubleDwi->SetPixel(index, pix); if (g==0) m_VolumeFractions.at(maxVolumeIndex+numFiberCompartments)->SetPixel(index, 1); } else { double extraAxonalVolume = m_VoxelVolume-intraAxonalVolume; // non-fiber volume if (extraAxonalVolume<0) { MITK_ERROR << "Coorupted intra-axonal signal voxel detected. Fiber volume larger voxel volume!"; extraAxonalVolume = 0; } double interAxonalVolume = 0; if (numFiberCompartments>1) interAxonalVolume = extraAxonalVolume * intraAxonalVolume/m_VoxelVolume; // inter-axonal fraction of non fiber compartment double other = extraAxonalVolume - interAxonalVolume; // rest of compartment if (other<0) { MITK_ERROR << "Coorupted signal voxel detected. Fiber volume larger voxel volume!"; other = 0; } // adjust non-fiber and intra-axonal signal for (int i=1; iGetPixel(index); if (intraAxonalVolume>0) // remove scaling by intra-axonal volume from inter-axonal compartment { if (g>=0) pix[g] /= intraAxonalVolume; else pix /= intraAxonalVolume; } if (m_Parameters.m_FiberModelList[i]->GetVolumeFractionImage()!=nullptr) { // ItkUcharImgType::IndexType maskIndex; // m_UpsampledMaskImage->TransformPhysicalPointToIndex(point, maskIndex); // if (!m_UpsampledMaskImage->GetLargestPossibleRegion().IsInside(maskIndex) || m_UpsampledMaskImage->GetPixel(maskIndex)==0) // continue; DoubleDwiType::IndexType newIndex; m_Parameters.m_FiberModelList[i]->GetVolumeFractionImage()->TransformPhysicalPointToIndex(point, newIndex); if (!m_Parameters.m_FiberModelList[i]->GetVolumeFractionImage()->GetLargestPossibleRegion().IsInside(newIndex)) continue; weight = m_Parameters.m_FiberModelList[i]->GetVolumeFractionImage()->GetPixel(newIndex)*m_VoxelVolume; } if (g>=0) pix[g] *= weight; else pix *= weight; m_CompartmentImages.at(i)->SetPixel(index, pix); if (g==0) m_VolumeFractions.at(i)->SetPixel(index, weight/m_VoxelVolume); } for (int i=0; iGetPixel(index); if (m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()!=nullptr) { // ItkUcharImgType::IndexType maskIndex; // m_UpsampledMaskImage->TransformPhysicalPointToIndex(point, maskIndex); // if (!m_UpsampledMaskImage->GetLargestPossibleRegion().IsInside(maskIndex) || m_UpsampledMaskImage->GetPixel(maskIndex)==0) // continue; DoubleDwiType::IndexType newIndex; m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->TransformPhysicalPointToIndex(point, newIndex); if (!m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetLargestPossibleRegion().IsInside(newIndex)) continue; weight = m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetPixel(newIndex)*m_VoxelVolume; } if (g>=0) pix[g] += m_Parameters.m_NonFiberModelList[i]->SimulateMeasurement(g)*weight; else pix += m_Parameters.m_NonFiberModelList[i]->SimulateMeasurement()*weight; m_CompartmentImages.at(i+numFiberCompartments)->SetPixel(index, pix); if (g==0) m_VolumeFractions.at(i+numFiberCompartments)->SetPixel(index, weight/m_VoxelVolume); } } } } template< class PixelType > itk::Point TractsToDWIImageFilter< PixelType >::GetItkPoint(double point[3]) { itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; return itkPoint; } template< class PixelType > itk::Vector TractsToDWIImageFilter< PixelType >::GetItkVector(double point[3]) { itk::Vector itkVector; itkVector[0] = point[0]; itkVector[1] = point[1]; itkVector[2] = point[2]; return itkVector; } template< class PixelType > vnl_vector_fixed TractsToDWIImageFilter< PixelType >::GetVnlVector(double point[3]) { vnl_vector_fixed vnlVector; vnlVector[0] = point[0]; vnlVector[1] = point[1]; vnlVector[2] = point[2]; return vnlVector; } template< class PixelType > vnl_vector_fixed TractsToDWIImageFilter< PixelType >::GetVnlVector(Vector& vector) { vnl_vector_fixed vnlVector; vnlVector[0] = vector[0]; vnlVector[1] = vector[1]; vnlVector[2] = vector[2]; return vnlVector; } template< class PixelType > double TractsToDWIImageFilter< PixelType >::RoundToNearest(double num) { return (num > 0.0) ? floor(num + 0.5) : ceil(num - 0.5); } template< class PixelType > std::string TractsToDWIImageFilter< PixelType >::GetTime() { m_TimeProbe.Stop(); unsigned long total = RoundToNearest(m_TimeProbe.GetTotal()); unsigned long hours = total/3600; unsigned long minutes = (total%3600)/60; unsigned long seconds = total%60; std::string out = ""; out.append(boost::lexical_cast(hours)); out.append(":"); out.append(boost::lexical_cast(minutes)); out.append(":"); out.append(boost::lexical_cast(seconds)); m_TimeProbe.Start(); return out; } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.h index 7e28f3403a..a04b540470 100755 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.h @@ -1,140 +1,144 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkTractsToDWIImageFilter_h__ #define __itkTractsToDWIImageFilter_h__ #include #include #include #include #include #include #include #include namespace itk { /** * \brief Generates artificial diffusion weighted image volume from the input fiberbundle using a generic multicompartment model. * See "Fiberfox: Facilitating the creation of realistic white matter software phantoms" (DOI: 10.1002/mrm.25045) for details. */ template< class PixelType > class TractsToDWIImageFilter : public ImageSource< itk::VectorImage< PixelType, 3 > > { public: typedef TractsToDWIImageFilter Self; typedef ImageSource< itk::VectorImage< PixelType, 3 > > Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; typedef typename Superclass::OutputImageType OutputImageType; typedef itk::Image ItkDoubleImgType; typedef itk::Image ItkFloatImgType; typedef itk::Image ItkUcharImgType; typedef mitk::FiberBundle::Pointer FiberBundleType; typedef itk::VectorImage< double, 3 > DoubleDwiType; typedef itk::Matrix MatrixType; typedef itk::Image< double, 2 > SliceType; typedef itk::VnlForwardFFTImageFilter::OutputImageType ComplexSliceType; typedef itk::Vector< double,3> DoubleVectorType; itkFactorylessNewMacro(Self) itkCloneMacro(Self) itkTypeMacro( TractsToDWIImageFilter, ImageSource ) /** Input */ itkSetMacro( FiberBundle, FiberBundleType ) ///< Input fiber bundle itkSetMacro( UseConstantRandSeed, bool ) ///< Seed for random generator. void SetParameters( FiberfoxParameters param ) ///< Simulation parameters. { m_Parameters = param; } /** Output */ FiberfoxParameters GetParameters(){ return m_Parameters; } std::vector< ItkDoubleImgType::Pointer > GetVolumeFractions() ///< one double image for each compartment containing the corresponding volume fraction per voxel { return m_VolumeFractions; } mitk::LevelWindow GetLevelWindow(){ return m_LevelWindow; } itkGetMacro( StatusText, std::string ) + itkGetMacro( PhaseImage, DoubleDwiType::Pointer ) + itkGetMacro( KspaceImage, DoubleDwiType::Pointer ) void GenerateData(); protected: TractsToDWIImageFilter(); virtual ~TractsToDWIImageFilter(); itk::Point GetItkPoint(double point[3]); itk::Vector GetItkVector(double point[3]); vnl_vector_fixed GetVnlVector(double point[3]); vnl_vector_fixed GetVnlVector(Vector< float, 3 >& vector); double RoundToNearest(double num); std::string GetTime(); /** Transform generated image compartment by compartment, channel by channel and slice by slice using DFT and add k-space artifacts. */ DoubleDwiType::Pointer DoKspaceStuff(std::vector< DoubleDwiType::Pointer >& images); /** Generate signal of non-fiber compartments. */ void SimulateExtraAxonalSignal(ItkUcharImgType::IndexType index, double intraAxonalVolume, int g=-1); /** Move fibers to simulate headmotion */ void SimulateMotion(int g=-1); void CheckVolumeFractionImages(); void InitializeData(); void InitializeFiberData(); // input mitk::FiberfoxParameters m_Parameters; FiberBundleType m_FiberBundle; // output typename OutputImageType::Pointer m_OutputImage; + typename DoubleDwiType::Pointer m_PhaseImage; + typename DoubleDwiType::Pointer m_KspaceImage; mitk::LevelWindow m_LevelWindow; std::vector< ItkDoubleImgType::Pointer > m_VolumeFractions; std::string m_StatusText; // MISC itk::TimeProbe m_TimeProbe; bool m_UseConstantRandSeed; bool m_MaskImageSet; ofstream m_Logfile; // signal generation FiberBundleType m_FiberBundleWorkingCopy; ///< we work on an upsampled version of the input bundle FiberBundleType m_FiberBundleTransformed; ///< transformed bundle simulating headmotion itk::Vector m_WorkingSpacing; itk::Point m_WorkingOrigin; ImageRegion<3> m_WorkingImageRegion; double m_VoxelVolume; std::vector< DoubleDwiType::Pointer > m_CompartmentImages; ItkUcharImgType::Pointer m_TransformedMaskImage; ///< copy of mask image (changes for each motion step) ItkUcharImgType::Pointer m_UpsampledMaskImage; ///< helper image for motion simulation DoubleVectorType m_Rotation; DoubleVectorType m_Translation; ImageRegion<3> m_CroppedRegion; double mmRadius; double segmentVolume; itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer m_RandGen; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkTractsToDWIImageFilter.cpp" #endif #endif diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h index 4f43a70a78..ec4b83c776 100644 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h @@ -1,325 +1,331 @@ /*=================================================================== 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 #include #include #include #include #include #include #include #include using namespace std; namespace mitk { /** Signal generation */ class SignalGenerationParameters { public: typedef itk::Image ItkDoubleImgType; typedef itk::Image ItkUcharImgType; typedef itk::Vector GradientType; typedef std::vector GradientListType; enum DiffusionDirectionMode { FIBER_TANGENT_DIRECTIONS, MAIN_FIBER_DIRECTIONS, RANDOM_DIRECTIONS }; SignalGenerationParameters() : m_SignalScale(100) , m_tEcho(100) , m_tRep(500) , m_tLine(1) , m_tInhom(50) , m_ReversePhase(false) + , m_PartialFourier(1.0) + , m_NoiseVariance(0.001) + , m_NumerOfCoils(1) , m_Bvalue(1000) , m_SimulateKspaceAcquisition(false) , m_AxonRadius(0) , m_DoDisablePartialVolume(false) , m_DiffusionDirectionMode(SignalGenerationParameters::FIBER_TANGENT_DIRECTIONS) , m_FiberSeparationThreshold(30) , m_Spikes(0) , m_SpikeAmplitude(1) , m_KspaceLineOffset(0) , m_EddyStrength(0) , m_Tau(70) , m_CroppingFactor(1) , m_DoAddGibbsRinging(false) , m_DoSimulateRelaxation(true) , m_DoAddMotion(false) , m_DoRandomizeMotion(true) , m_FrequencyMap(NULL) , m_MaskImage(NULL) { m_ImageRegion.SetSize(0, 12); m_ImageRegion.SetSize(1, 12); m_ImageRegion.SetSize(2, 3); m_ImageSpacing.Fill(2.0); m_ImageOrigin.Fill(0.0); m_ImageDirection.SetIdentity(); m_Translation.Fill(0.0); m_Rotation.Fill(0.0); SetNumWeightedVolumes(6); } /** input/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 */ double m_SignalScale; ///< Scaling factor for output signal (before noise is added). double m_tEcho; ///< Echo time TE. double m_tRep; ///< Echo time TR. double m_tLine; ///< k-space line readout time (dwell time). double m_tInhom; ///< T2' bool m_ReversePhase; ///< If true, the phase readout direction will be inverted (-y instead of y) + double m_PartialFourier; ///< Partial fourier factor (0.5-1) + double m_NoiseVariance; + unsigned int m_NumerOfCoils; double m_Bvalue; ///< Acquisition b-value bool m_SimulateKspaceAcquisition;///< Flag to enable/disable k-space acquisition simulation double m_AxonRadius; ///< Determines compartment volume fractions (0 == automatic axon radius estimation) bool m_DoDisablePartialVolume; ///< Disable partial volume effects. Each voxel is either all fiber or all non-fiber. DiffusionDirectionMode m_DiffusionDirectionMode; ///< Determines how the main diffusion direction of the signal models is selected double m_FiberSeparationThreshold; ///< Used for random and and mein fiber deriction DiffusionDirectionMode /** Artifacts and other effects */ unsigned 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 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_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. 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. inline void GenerateGradientHalfShell(); ///< Generates half shell of gradient directions (with m_NumGradients non-zero directions) inline std::vector< int > GetBaselineIndices(); ///< Returns list of nun-diffusion-weighted image volume indices inline unsigned int GetFirstBaselineIndex(); ///< Returns index of first non-diffusion-weighted image volume inline bool IsBaselineIndex(unsigned int idx); ///< Checks if image volume with given index is non-diffusion-weighted volume or not. inline unsigned int GetNumWeightedVolumes(); ///< Get number of diffusion-weighted image volumes inline unsigned int GetNumBaselineVolumes(); ///< Get number of non-diffusion-weighted image volumes inline unsigned int GetNumVolumes(); ///< Get number of baseline and diffusion-weighted image volumes inline GradientListType GetGradientDirections(); ///< Return gradient direction container inline GradientType GetGradientDirection(unsigned int i); inline void SetNumWeightedVolumes(int numGradients); ///< Automaticall calls GenerateGradientHalfShell() afterwards. inline void SetGradienDirections(GradientListType gradientList); inline void SetGradienDirections(mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradientList); protected: unsigned int m_NumGradients; ///< Number of diffusion-weighted image volumes. unsigned int m_NumBaseline; ///< Number of non-diffusion-weighted image volumes. GradientListType m_GradientDirections; ///< Total number of image volumes. }; /** Fiber generation */ class FiberGenerationParameters { public: enum FiberDistribution{ DISTRIBUTE_UNIFORM, // distribute fibers uniformly in the ROIs DISTRIBUTE_GAUSSIAN // distribute fibers using a 2D gaussian }; typedef vector< vector< mitk::PlanarEllipse::Pointer > > FiducialListType; typedef vector< vector< unsigned int > > FlipListType; FiberGenerationParameters() : m_Distribution(DISTRIBUTE_UNIFORM) , m_Density(100) , m_Variance(100) , m_Sampling(1) , m_Tension(0) , m_Continuity(0) , m_Bias(0) { m_Rotation.Fill(0.0); m_Translation.Fill(0.0); m_Scale.Fill(1.0); } FiberDistribution m_Distribution; unsigned int m_Density; double m_Variance; double m_Sampling; double m_Tension; double m_Continuity; double m_Bias; mitk::Vector3D m_Rotation; mitk::Vector3D m_Translation; mitk::Vector3D m_Scale; FlipListType m_FlipList; ///< contains flags indicating a flip of the 2D fiber x-coordinates (needed to resolve some unwanted fiber twisting) FiducialListType m_Fiducials; ///< container of the planar ellipses used as fiducials for the fiber generation process }; /** GUI persistence, input, output, ... */ class MiscFiberfoxParameters { public: MiscFiberfoxParameters() : m_ResultNode(DataNode::New()) , m_ParentNode(NULL) , m_SignalModelString("") , m_ArtifactModelString("") , m_OutputPath("") , m_CheckOutputVolumeFractionsBox(false) , m_CheckAdvancedSignalOptionsBox(false) , m_CheckAddNoiseBox(false) , m_CheckAddGhostsBox(false) , m_CheckAddAliasingBox(false) , m_CheckAddSpikesBox(false) , m_CheckAddEddyCurrentsBox(false) , m_CheckAddDistortionsBox(false) , m_CheckRealTimeFibersBox(true) , m_CheckAdvancedFiberOptionsBox(false) , m_CheckConstantRadiusBox(false) , m_CheckIncludeFiducialsBox(true) {} DataNode::Pointer m_ResultNode; ///< Stores resulting image. DataNode::Pointer m_ParentNode; ///< Parent node of 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. /** member variables that store the check-state of GUI checkboxes */ // image generation bool m_CheckOutputVolumeFractionsBox; bool m_CheckAdvancedSignalOptionsBox; bool m_CheckAddNoiseBox; bool m_CheckAddGhostsBox; bool m_CheckAddAliasingBox; bool m_CheckAddSpikesBox; bool m_CheckAddEddyCurrentsBox; bool m_CheckAddDistortionsBox; // fiber generation bool m_CheckRealTimeFibersBox; bool m_CheckAdvancedFiberOptionsBox; bool m_CheckConstantRadiusBox; bool m_CheckIncludeFiducialsBox; }; /** * \brief Datastructure to manage the Fiberfox signal generation parameters. * */ template< class ScalarType = double > class FiberfoxParameters { public: typedef itk::Image ItkDoubleImgType; typedef itk::Image ItkUcharImgType; typedef DiffusionSignalModel DiffusionModelType; typedef std::vector< DiffusionModelType* > DiffusionModelListType; typedef DiffusionNoiseModel NoiseModelType; FiberfoxParameters(); ~FiberfoxParameters(); /** Get same parameter object with different template parameter */ template< class OutType > FiberfoxParameters< OutType > CopyParameters() { FiberfoxParameters< OutType > out; out.m_FiberGen = m_FiberGen; out.m_SignalGen = m_SignalGen; out.m_Misc = m_Misc; 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()); } for (unsigned int i=0; i* outModel = NULL; mitk::DiffusionSignalModel* signalModel = NULL; if (i*>(signalModel)) outModel = new mitk::StickModel(dynamic_cast*>(signalModel)); else if (dynamic_cast*>(signalModel)) outModel = new mitk::TensorModel(dynamic_cast*>(signalModel)); else if (dynamic_cast*>(signalModel)) outModel = new mitk::RawShModel(dynamic_cast*>(signalModel)); else if (dynamic_cast*>(signalModel)) outModel = new mitk::BallModel(dynamic_cast*>(signalModel)); else if (dynamic_cast*>(signalModel)) outModel = new mitk::AstroStickModel(dynamic_cast*>(signalModel)); else if (dynamic_cast*>(signalModel)) outModel = new mitk::DotModel(dynamic_cast*>(signalModel)); if (i // Blueberry #include #include // Qmitk #include "QmitkFiberfoxView.h" // MITK #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "usModuleRegistry.h" #include #include #include #include #include #include #include #include #include #include "mitkNodePredicateDataType.h" #include #include #include #include #define _USE_MATH_DEFINES #include QmitkFiberfoxWorker::QmitkFiberfoxWorker(QmitkFiberfoxView* view) : m_View(view) { } void QmitkFiberfoxWorker::run() { try{ switch (m_FilterType) { case 0: m_View->m_TractsToDwiFilter->Update(); break; case 1: m_View->m_ArtifactsToDwiFilter->Update(); break; } } catch( ... ) { } m_View->m_Thread.quit(); } const std::string QmitkFiberfoxView::VIEW_ID = "org.mitk.views.fiberfoxview"; QmitkFiberfoxView::QmitkFiberfoxView() : QmitkAbstractView() , m_Controls( 0 ) , m_SelectedImageNode( NULL ) , m_Worker(this) , m_ThreadIsRunning(false) { m_Worker.moveToThread(&m_Thread); connect(&m_Thread, SIGNAL(started()), this, SLOT(BeforeThread())); connect(&m_Thread, SIGNAL(started()), &m_Worker, SLOT(run())); connect(&m_Thread, SIGNAL(finished()), this, SLOT(AfterThread())); connect(&m_Thread, SIGNAL(terminated()), this, SLOT(AfterThread())); m_SimulationTimer = new QTimer(this); } void QmitkFiberfoxView::KillThread() { MITK_INFO << "Aborting DWI simulation."; switch (m_Worker.m_FilterType) { case 0: m_TractsToDwiFilter->SetAbortGenerateData(true); break; case 1: m_ArtifactsToDwiFilter->SetAbortGenerateData(true); break; } m_Controls->m_AbortSimulationButton->setEnabled(false); m_Controls->m_AbortSimulationButton->setText("Aborting simulation ..."); } void QmitkFiberfoxView::BeforeThread() { m_SimulationTime = QTime::currentTime(); m_SimulationTimer->start(100); m_Controls->m_AbortSimulationButton->setVisible(true); m_Controls->m_GenerateImageButton->setVisible(false); m_Controls->m_SimulationStatusText->setVisible(true); m_ThreadIsRunning = true; } void QmitkFiberfoxView::AfterThread() { UpdateSimulationStatus(); m_SimulationTimer->stop(); m_Controls->m_AbortSimulationButton->setVisible(false); m_Controls->m_AbortSimulationButton->setEnabled(true); m_Controls->m_AbortSimulationButton->setText("Abort simulation"); m_Controls->m_GenerateImageButton->setVisible(true); m_ThreadIsRunning = false; QString statusText; FiberfoxParameters parameters; mitk::Image::Pointer mitkImage = mitk::Image::New(); switch (m_Worker.m_FilterType) { case 0: { statusText = QString(m_TractsToDwiFilter->GetStatusText().c_str()); if (m_TractsToDwiFilter->GetAbortGenerateData()) { MITK_INFO << "Simulation aborted."; return; } parameters = m_TractsToDwiFilter->GetParameters(); mitkImage = mitk::GrabItkImageMemory( m_TractsToDwiFilter->GetOutput() ); mitkImage->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( parameters.m_SignalGen.GetGradientDirections() )); mitkImage->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( parameters.m_SignalGen.m_Bvalue )); mitk::DiffusionPropertyHelper propertyHelper( mitkImage ); propertyHelper.InitializeImage(); parameters.m_Misc.m_ResultNode->SetData( mitkImage ); parameters.m_Misc.m_ResultNode->SetName(parameters.m_Misc.m_ParentNode->GetName() +"_D"+QString::number(parameters.m_SignalGen.m_ImageRegion.GetSize(0)).toStdString() +"-"+QString::number(parameters.m_SignalGen.m_ImageRegion.GetSize(1)).toStdString() +"-"+QString::number(parameters.m_SignalGen.m_ImageRegion.GetSize(2)).toStdString() +"_S"+QString::number(parameters.m_SignalGen.m_ImageSpacing[0]).toStdString() +"-"+QString::number(parameters.m_SignalGen.m_ImageSpacing[1]).toStdString() +"-"+QString::number(parameters.m_SignalGen.m_ImageSpacing[2]).toStdString() +"_b"+QString::number(parameters.m_SignalGen.m_Bvalue).toStdString() +"_"+parameters.m_Misc.m_SignalModelString +parameters.m_Misc.m_ArtifactModelString); GetDataStorage()->Add(parameters.m_Misc.m_ResultNode, parameters.m_Misc.m_ParentNode); parameters.m_Misc.m_ResultNode->SetProperty( "levelwindow", mitk::LevelWindowProperty::New(m_TractsToDwiFilter->GetLevelWindow()) ); if (m_Controls->m_VolumeFractionsBox->isChecked()) { std::vector< itk::TractsToDWIImageFilter< short >::ItkDoubleImgType::Pointer > volumeFractions = m_TractsToDwiFilter->GetVolumeFractions(); for (unsigned int k=0; kInitializeByItk(volumeFractions.at(k).GetPointer()); image->SetVolume(volumeFractions.at(k)->GetBufferPointer()); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); - node->SetName(parameters.m_Misc.m_ParentNode->GetName()+"_CompartmentVolume-"+QString::number(k).toStdString()); - GetDataStorage()->Add(node, parameters.m_Misc.m_ParentNode); + node->SetName("CompartmentVolume-"+QString::number(k).toStdString()); + GetDataStorage()->Add(node, parameters.m_Misc.m_ResultNode); + } + + if (m_TractsToDwiFilter->GetPhaseImage().IsNotNull()) + { + mitk::Image::Pointer phaseImage = mitk::Image::New(); + itk::TractsToDWIImageFilter< short >::DoubleDwiType::Pointer itkPhase = m_TractsToDwiFilter->GetPhaseImage(); + phaseImage = mitk::GrabItkImageMemory( itkPhase.GetPointer() ); + mitk::DataNode::Pointer phaseNode = mitk::DataNode::New(); + phaseNode->SetData( phaseImage ); + phaseNode->SetName("Phase Image"); + GetDataStorage()->Add(phaseNode, parameters.m_Misc.m_ResultNode); + } + + if (m_TractsToDwiFilter->GetKspaceImage().IsNotNull()) + { + mitk::Image::Pointer image = mitk::Image::New(); + itk::TractsToDWIImageFilter< short >::DoubleDwiType::Pointer itkImage = m_TractsToDwiFilter->GetKspaceImage(); + image = mitk::GrabItkImageMemory( itkImage.GetPointer() ); + mitk::DataNode::Pointer node = mitk::DataNode::New(); + node->SetData( image ); + node->SetName("k-Space"); + GetDataStorage()->Add(node, parameters.m_Misc.m_ResultNode); } } m_TractsToDwiFilter = NULL; break; } case 1: { statusText = QString(m_ArtifactsToDwiFilter->GetStatusText().c_str()); if (m_ArtifactsToDwiFilter->GetAbortGenerateData()) { MITK_INFO << "Simulation aborted."; return; } parameters = m_ArtifactsToDwiFilter->GetParameters().CopyParameters(); mitk::Image::Pointer diffImg = dynamic_cast(parameters.m_Misc.m_ParentNode->GetData()); mitkImage = mitk::GrabItkImageMemory( m_ArtifactsToDwiFilter->GetOutput() ); mitkImage->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( static_cast( diffImg->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer() ) ); mitkImage->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( static_cast( diffImg->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue() ) ); mitk::DiffusionPropertyHelper propertyHelper( mitkImage ); propertyHelper.InitializeImage(); parameters.m_Misc.m_ResultNode->SetData( mitkImage ); parameters.m_Misc.m_ResultNode->SetName(parameters.m_Misc.m_ParentNode->GetName()+parameters.m_Misc.m_ArtifactModelString); GetDataStorage()->Add(parameters.m_Misc.m_ResultNode, parameters.m_Misc.m_ParentNode); m_ArtifactsToDwiFilter = NULL; break; } } mitk::BaseData::Pointer basedata = parameters.m_Misc.m_ResultNode->GetData(); if (basedata.IsNotNull()) { mitk::RenderingManager::GetInstance()->InitializeViews( basedata->GetTimeGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true ); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } if (!parameters.m_Misc.m_OutputPath.empty()) { try{ QString outputFileName(parameters.m_Misc.m_OutputPath.c_str()); outputFileName += parameters.m_Misc.m_ResultNode->GetName().c_str(); outputFileName.replace(QString("."), QString("_")); outputFileName += ".dwi"; QString status("Saving output image to "); status += outputFileName; m_Controls->m_SimulationStatusText->append(status); mitk::IOUtil::SaveBaseData(mitkImage, outputFileName.toStdString()); m_Controls->m_SimulationStatusText->append("File saved successfully."); } catch (itk::ExceptionObject &e) { QString status("Exception during DWI writing: "); status += e.GetDescription(); m_Controls->m_SimulationStatusText->append(status); } catch (...) { m_Controls->m_SimulationStatusText->append("Unknown exception during DWI writing!"); } } parameters.m_SignalGen.m_FrequencyMap = NULL; } void QmitkFiberfoxView::UpdateSimulationStatus() { QString statusText; switch (m_Worker.m_FilterType) { case 0: statusText = QString(m_TractsToDwiFilter->GetStatusText().c_str()); break; case 1: statusText = QString(m_ArtifactsToDwiFilter->GetStatusText().c_str()); break; } if (QString::compare(m_SimulationStatusText,statusText)!=0) { m_Controls->m_SimulationStatusText->clear(); statusText = "
"+statusText+"
"; m_Controls->m_SimulationStatusText->setText(statusText); QScrollBar *vScrollBar = m_Controls->m_SimulationStatusText->verticalScrollBar(); vScrollBar->triggerAction(QScrollBar::SliderToMaximum); } } // Destructor QmitkFiberfoxView::~QmitkFiberfoxView() { delete m_SimulationTimer; } void QmitkFiberfoxView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkFiberfoxViewControls; m_Controls->setupUi( parent ); - // commented out - m_Controls->m_DiffusionDirectionBox->setVisible(false); - m_Controls->label_3->setVisible(false); - m_Controls->m_SeparationAngleBox->setVisible(false); - m_Controls->label_4->setVisible(false); - // - m_Controls->m_StickWidget1->setVisible(true); m_Controls->m_StickWidget2->setVisible(false); m_Controls->m_ZeppelinWidget1->setVisible(false); m_Controls->m_ZeppelinWidget2->setVisible(false); m_Controls->m_TensorWidget1->setVisible(false); m_Controls->m_TensorWidget2->setVisible(false); m_Controls->m_BallWidget1->setVisible(true); m_Controls->m_BallWidget2->setVisible(false); m_Controls->m_AstrosticksWidget1->setVisible(false); m_Controls->m_AstrosticksWidget2->setVisible(false); m_Controls->m_DotWidget1->setVisible(false); m_Controls->m_DotWidget2->setVisible(false); m_Controls->m_PrototypeWidget1->setVisible(false); m_Controls->m_PrototypeWidget2->setVisible(false); m_Controls->m_PrototypeWidget3->setVisible(false); m_Controls->m_PrototypeWidget4->setVisible(false); m_Controls->m_PrototypeWidget3->SetMinFa(0.0); m_Controls->m_PrototypeWidget3->SetMaxFa(0.15); m_Controls->m_PrototypeWidget4->SetMinFa(0.0); m_Controls->m_PrototypeWidget4->SetMaxFa(0.15); m_Controls->m_PrototypeWidget3->SetMinAdc(0.0); m_Controls->m_PrototypeWidget3->SetMaxAdc(0.001); m_Controls->m_PrototypeWidget4->SetMinAdc(0.003); m_Controls->m_PrototypeWidget4->SetMaxAdc(0.004); m_Controls->m_Comp2FractionFrame->setVisible(false); m_Controls->m_Comp4FractionFrame->setVisible(false); m_Controls->m_DiffusionPropsMessage->setVisible(false); m_Controls->m_GeometryMessage->setVisible(false); m_Controls->m_AdvancedSignalOptionsFrame->setVisible(false); m_Controls->m_AdvancedFiberOptionsFrame->setVisible(false); m_Controls->m_VarianceBox->setVisible(false); m_Controls->m_NoiseFrame->setVisible(false); m_Controls->m_GhostFrame->setVisible(false); m_Controls->m_DistortionsFrame->setVisible(false); m_Controls->m_EddyFrame->setVisible(false); m_Controls->m_SpikeFrame->setVisible(false); m_Controls->m_AliasingFrame->setVisible(false); m_Controls->m_MotionArtifactFrame->setVisible(false); m_ParameterFile = QDir::currentPath()+"/param.ffp"; m_Controls->m_AbortSimulationButton->setVisible(false); m_Controls->m_SimulationStatusText->setVisible(false); m_Controls->m_FrequencyMapBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_Comp1VolumeFraction->SetDataStorage(this->GetDataStorage()); m_Controls->m_Comp2VolumeFraction->SetDataStorage(this->GetDataStorage()); m_Controls->m_Comp3VolumeFraction->SetDataStorage(this->GetDataStorage()); m_Controls->m_Comp4VolumeFraction->SetDataStorage(this->GetDataStorage()); m_Controls->m_MaskComboBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_TemplateComboBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_FiberBundleComboBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isFiberBundle = mitk::TNodePredicateDataType::New(); mitk::TNodePredicateDataType::Pointer isMitkImage = mitk::TNodePredicateDataType::New(); mitk::NodePredicateIsDWI::Pointer isDwi = mitk::NodePredicateIsDWI::New( ); mitk::NodePredicateDataType::Pointer isDti = mitk::NodePredicateDataType::New("TensorImage"); mitk::NodePredicateDataType::Pointer isQbi = mitk::NodePredicateDataType::New("QBallImage"); mitk::NodePredicateOr::Pointer isDiffusionImage = mitk::NodePredicateOr::New(isDwi, isDti); isDiffusionImage = mitk::NodePredicateOr::New(isDiffusionImage, isQbi); mitk::NodePredicateNot::Pointer noDiffusionImage = mitk::NodePredicateNot::New(isDiffusionImage); mitk::NodePredicateAnd::Pointer isNonDiffMitkImage = mitk::NodePredicateAnd::New(isMitkImage, noDiffusionImage); mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); mitk::NodePredicateAnd::Pointer isBinaryMitkImage = mitk::NodePredicateAnd::New( isNonDiffMitkImage, isBinaryPredicate ); m_Controls->m_FrequencyMapBox->SetPredicate(isNonDiffMitkImage); m_Controls->m_Comp1VolumeFraction->SetPredicate(isNonDiffMitkImage); m_Controls->m_Comp1VolumeFraction->SetZeroEntryText("--"); m_Controls->m_Comp2VolumeFraction->SetPredicate(isNonDiffMitkImage); m_Controls->m_Comp2VolumeFraction->SetZeroEntryText("--"); m_Controls->m_Comp3VolumeFraction->SetPredicate(isNonDiffMitkImage); m_Controls->m_Comp3VolumeFraction->SetZeroEntryText("--"); m_Controls->m_Comp4VolumeFraction->SetPredicate(isNonDiffMitkImage); m_Controls->m_Comp4VolumeFraction->SetZeroEntryText("--"); m_Controls->m_MaskComboBox->SetPredicate(isBinaryMitkImage); m_Controls->m_MaskComboBox->SetZeroEntryText("--"); m_Controls->m_TemplateComboBox->SetPredicate(isMitkImage); m_Controls->m_TemplateComboBox->SetZeroEntryText("--"); m_Controls->m_FiberBundleComboBox->SetPredicate(isFiberBundle); m_Controls->m_FiberBundleComboBox->SetZeroEntryText("--"); // mitk::NodePredicateDimension::Pointer dimensionPredicate = mitk::NodePredicateDimension::New(3); connect( m_SimulationTimer, SIGNAL(timeout()), this, SLOT(UpdateSimulationStatus()) ); connect((QObject*) m_Controls->m_AbortSimulationButton, SIGNAL(clicked()), (QObject*) this, SLOT(KillThread())); connect((QObject*) m_Controls->m_GenerateImageButton, SIGNAL(clicked()), (QObject*) this, SLOT(GenerateImage())); connect((QObject*) m_Controls->m_GenerateFibersButton, SIGNAL(clicked()), (QObject*) this, SLOT(GenerateFibers())); connect((QObject*) m_Controls->m_CircleButton, SIGNAL(clicked()), (QObject*) this, SLOT(OnDrawROI())); connect((QObject*) m_Controls->m_FlipButton, SIGNAL(clicked()), (QObject*) this, SLOT(OnFlipButton())); connect((QObject*) m_Controls->m_JoinBundlesButton, SIGNAL(clicked()), (QObject*) this, SLOT(JoinBundles())); connect((QObject*) m_Controls->m_VarianceBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnVarianceChanged(double))); connect((QObject*) m_Controls->m_DistributionBox, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(OnDistributionChanged(int))); connect((QObject*) m_Controls->m_FiberDensityBox, SIGNAL(valueChanged(int)), (QObject*) this, SLOT(OnFiberDensityChanged(int))); connect((QObject*) m_Controls->m_FiberSamplingBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnFiberSamplingChanged(double))); connect((QObject*) m_Controls->m_TensionBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnTensionChanged(double))); connect((QObject*) m_Controls->m_ContinuityBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnContinuityChanged(double))); connect((QObject*) m_Controls->m_BiasBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnBiasChanged(double))); connect((QObject*) m_Controls->m_AddNoise, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddNoise(int))); connect((QObject*) m_Controls->m_AddGhosts, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddGhosts(int))); connect((QObject*) m_Controls->m_AddDistortions, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddDistortions(int))); connect((QObject*) m_Controls->m_AddEddy, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddEddy(int))); connect((QObject*) m_Controls->m_AddSpikes, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddSpikes(int))); connect((QObject*) m_Controls->m_AddAliasing, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddAliasing(int))); connect((QObject*) m_Controls->m_AddMotion, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddMotion(int))); connect((QObject*) m_Controls->m_ConstantRadiusBox, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnConstantRadius(int))); connect((QObject*) m_Controls->m_CopyBundlesButton, SIGNAL(clicked()), (QObject*) this, SLOT(CopyBundles())); connect((QObject*) m_Controls->m_TransformBundlesButton, SIGNAL(clicked()), (QObject*) this, SLOT(ApplyTransform())); connect((QObject*) m_Controls->m_AlignOnGrid, SIGNAL(clicked()), (QObject*) this, SLOT(AlignOnGrid())); connect((QObject*) m_Controls->m_Compartment1Box, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(Comp1ModelFrameVisibility(int))); connect((QObject*) m_Controls->m_Compartment2Box, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(Comp2ModelFrameVisibility(int))); connect((QObject*) m_Controls->m_Compartment3Box, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(Comp3ModelFrameVisibility(int))); connect((QObject*) m_Controls->m_Compartment4Box, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(Comp4ModelFrameVisibility(int))); connect((QObject*) m_Controls->m_AdvancedOptionsBox, SIGNAL( stateChanged(int)), (QObject*) this, SLOT(ShowAdvancedOptions(int))); connect((QObject*) m_Controls->m_AdvancedOptionsBox_2, SIGNAL( stateChanged(int)), (QObject*) this, SLOT(ShowAdvancedOptions(int))); connect((QObject*) m_Controls->m_SaveParametersButton, SIGNAL(clicked()), (QObject*) this, SLOT(SaveParameters())); connect((QObject*) m_Controls->m_LoadParametersButton, SIGNAL(clicked()), (QObject*) this, SLOT(LoadParameters())); connect((QObject*) m_Controls->m_OutputPathButton, SIGNAL(clicked()), (QObject*) this, SLOT(SetOutputPath())); connect((QObject*) m_Controls->m_MaskComboBox, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(OnMaskSelected(int))); connect((QObject*) m_Controls->m_TemplateComboBox, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(OnTemplateSelected(int))); connect((QObject*) m_Controls->m_FiberBundleComboBox, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(OnFibSelected(int))); } } void QmitkFiberfoxView::OnMaskSelected(int value) { UpdateGui(); } void QmitkFiberfoxView::OnTemplateSelected(int value) { UpdateGui(); } void QmitkFiberfoxView::OnFibSelected(int value) { UpdateGui(); } template< class ScalarType > FiberfoxParameters< ScalarType > QmitkFiberfoxView::UpdateImageParameters() { FiberfoxParameters< ScalarType > parameters; parameters.m_Misc.m_OutputPath = ""; parameters.m_Misc.m_CheckAdvancedFiberOptionsBox = m_Controls->m_AdvancedOptionsBox->isChecked(); parameters.m_Misc.m_CheckAdvancedSignalOptionsBox = m_Controls->m_AdvancedOptionsBox_2->isChecked(); parameters.m_Misc.m_CheckOutputVolumeFractionsBox = m_Controls->m_VolumeFractionsBox->isChecked(); string outputPath = m_Controls->m_SavePathEdit->text().toStdString(); if (outputPath.compare("-")!=0) { parameters.m_Misc.m_OutputPath = outputPath; parameters.m_Misc.m_OutputPath += "/"; } if (m_Controls->m_MaskComboBox->GetSelectedNode().IsNotNull()) { mitk::Image::Pointer mitkMaskImage = dynamic_cast(m_Controls->m_MaskComboBox->GetSelectedNode()->GetData()); mitk::CastToItkImage(mitkMaskImage, parameters.m_SignalGen.m_MaskImage); itk::ImageDuplicator::Pointer duplicator = itk::ImageDuplicator::New(); duplicator->SetInputImage(parameters.m_SignalGen.m_MaskImage); duplicator->Update(); parameters.m_SignalGen.m_MaskImage = duplicator->GetOutput(); } if (m_Controls->m_TemplateComboBox->GetSelectedNode().IsNotNull() && mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( m_Controls->m_TemplateComboBox->GetSelectedNode())) // use parameters of selected DWI { mitk::Image::Pointer dwi = dynamic_cast(m_Controls->m_TemplateComboBox->GetSelectedNode()->GetData()); ItkDwiType::Pointer itkVectorImagePointer = ItkDwiType::New(); mitk::CastToItkImage(dwi, itkVectorImagePointer); parameters.m_SignalGen.m_ImageRegion = itkVectorImagePointer->GetLargestPossibleRegion(); parameters.m_SignalGen.m_ImageSpacing = itkVectorImagePointer->GetSpacing(); parameters.m_SignalGen.m_ImageOrigin = itkVectorImagePointer->GetOrigin(); parameters.m_SignalGen.m_ImageDirection = itkVectorImagePointer->GetDirection(); parameters.m_SignalGen.m_Bvalue = static_cast(dwi->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue(); parameters.m_SignalGen.SetGradienDirections(static_cast( dwi->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer()); } else if (m_Controls->m_TemplateComboBox->GetSelectedNode().IsNotNull()) // use geometry of selected image { mitk::Image::Pointer img = dynamic_cast(m_Controls->m_TemplateComboBox->GetSelectedNode()->GetData()); itk::Image< float, 3 >::Pointer itkImg = itk::Image< float, 3 >::New(); CastToItkImage< itk::Image< float, 3 > >(img, itkImg); parameters.m_SignalGen.m_ImageRegion = itkImg->GetLargestPossibleRegion(); parameters.m_SignalGen.m_ImageSpacing = itkImg->GetSpacing(); parameters.m_SignalGen.m_ImageOrigin = itkImg->GetOrigin(); parameters.m_SignalGen.m_ImageDirection = itkImg->GetDirection(); parameters.m_SignalGen.SetNumWeightedVolumes(m_Controls->m_NumGradientsBox->value()); parameters.m_SignalGen.m_Bvalue = m_Controls->m_BvalueBox->value(); } else if (parameters.m_SignalGen.m_MaskImage.IsNotNull()) // use geometry of mask image { ItkUcharImgType::Pointer itkImg = parameters.m_SignalGen.m_MaskImage; parameters.m_SignalGen.m_ImageRegion = itkImg->GetLargestPossibleRegion(); parameters.m_SignalGen.m_ImageSpacing = itkImg->GetSpacing(); parameters.m_SignalGen.m_ImageOrigin = itkImg->GetOrigin(); parameters.m_SignalGen.m_ImageDirection = itkImg->GetDirection(); parameters.m_SignalGen.SetNumWeightedVolumes(m_Controls->m_NumGradientsBox->value()); parameters.m_SignalGen.m_Bvalue = m_Controls->m_BvalueBox->value(); } else // use GUI parameters { parameters.m_SignalGen.m_ImageRegion.SetSize(0, m_Controls->m_SizeX->value()); parameters.m_SignalGen.m_ImageRegion.SetSize(1, m_Controls->m_SizeY->value()); parameters.m_SignalGen.m_ImageRegion.SetSize(2, m_Controls->m_SizeZ->value()); parameters.m_SignalGen.m_ImageSpacing[0] = m_Controls->m_SpacingX->value(); parameters.m_SignalGen.m_ImageSpacing[1] = m_Controls->m_SpacingY->value(); parameters.m_SignalGen.m_ImageSpacing[2] = m_Controls->m_SpacingZ->value(); parameters.m_SignalGen.m_ImageOrigin[0] = parameters.m_SignalGen.m_ImageSpacing[0]/2; parameters.m_SignalGen.m_ImageOrigin[1] = parameters.m_SignalGen.m_ImageSpacing[1]/2; parameters.m_SignalGen.m_ImageOrigin[2] = parameters.m_SignalGen.m_ImageSpacing[2]/2; parameters.m_SignalGen.m_ImageDirection.SetIdentity(); parameters.m_SignalGen.SetNumWeightedVolumes(m_Controls->m_NumGradientsBox->value()); parameters.m_SignalGen.m_Bvalue = m_Controls->m_BvalueBox->value(); parameters.m_SignalGen.GenerateGradientHalfShell(); } // signal relaxation parameters.m_SignalGen.m_DoSimulateRelaxation = m_Controls->m_RelaxationBox->isChecked(); parameters.m_SignalGen.m_SimulateKspaceAcquisition = parameters.m_SignalGen.m_DoSimulateRelaxation; if (parameters.m_SignalGen.m_DoSimulateRelaxation && m_Controls->m_FiberBundleComboBox->GetSelectedNode().IsNotNull() ) parameters.m_Misc.m_ArtifactModelString += "_RELAX"; // N/2 ghosts parameters.m_Misc.m_CheckAddGhostsBox = m_Controls->m_AddGhosts->isChecked(); if (m_Controls->m_AddGhosts->isChecked()) { parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; parameters.m_Misc.m_ArtifactModelString += "_GHOST"; parameters.m_SignalGen.m_KspaceLineOffset = m_Controls->m_kOffsetBox->value(); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Ghost", DoubleProperty::New(parameters.m_SignalGen.m_KspaceLineOffset)); } else parameters.m_SignalGen.m_KspaceLineOffset = 0; // Aliasing parameters.m_Misc.m_CheckAddAliasingBox = m_Controls->m_AddAliasing->isChecked(); if (m_Controls->m_AddAliasing->isChecked()) { parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; parameters.m_Misc.m_ArtifactModelString += "_ALIASING"; parameters.m_SignalGen.m_CroppingFactor = (100-m_Controls->m_WrapBox->value())/100; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Aliasing", DoubleProperty::New(m_Controls->m_WrapBox->value())); } // Spikes parameters.m_Misc.m_CheckAddSpikesBox = m_Controls->m_AddSpikes->isChecked(); if (m_Controls->m_AddSpikes->isChecked()) { parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; parameters.m_SignalGen.m_Spikes = m_Controls->m_SpikeNumBox->value(); parameters.m_SignalGen.m_SpikeAmplitude = m_Controls->m_SpikeScaleBox->value(); parameters.m_Misc.m_ArtifactModelString += "_SPIKES"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Spikes.Number", IntProperty::New(parameters.m_SignalGen.m_Spikes)); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Spikes.Amplitude", DoubleProperty::New(parameters.m_SignalGen.m_SpikeAmplitude)); } // gibbs ringing parameters.m_SignalGen.m_DoAddGibbsRinging = m_Controls->m_AddGibbsRinging->isChecked(); if (m_Controls->m_AddGibbsRinging->isChecked()) { parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Ringing", BoolProperty::New(true)); parameters.m_Misc.m_ArtifactModelString += "_RINGING"; } // add distortions parameters.m_Misc.m_CheckAddDistortionsBox = m_Controls->m_AddDistortions->isChecked(); if (m_Controls->m_AddDistortions->isChecked() && m_Controls->m_FrequencyMapBox->GetSelectedNode().IsNotNull()) { mitk::DataNode::Pointer fMapNode = m_Controls->m_FrequencyMapBox->GetSelectedNode(); mitk::Image* img = dynamic_cast(fMapNode->GetData()); ItkDoubleImgType::Pointer itkImg = ItkDoubleImgType::New(); CastToItkImage< ItkDoubleImgType >(img, itkImg); if (m_Controls->m_TemplateComboBox->GetSelectedNode().IsNull()) // use geometry of frequency map { parameters.m_SignalGen.m_ImageRegion = itkImg->GetLargestPossibleRegion(); parameters.m_SignalGen.m_ImageSpacing = itkImg->GetSpacing(); parameters.m_SignalGen.m_ImageOrigin = itkImg->GetOrigin(); parameters.m_SignalGen.m_ImageDirection = itkImg->GetDirection(); } if (parameters.m_SignalGen.m_ImageRegion.GetSize(0)==itkImg->GetLargestPossibleRegion().GetSize(0) && parameters.m_SignalGen.m_ImageRegion.GetSize(1)==itkImg->GetLargestPossibleRegion().GetSize(1) && parameters.m_SignalGen.m_ImageRegion.GetSize(2)==itkImg->GetLargestPossibleRegion().GetSize(2)) { parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; itk::ImageDuplicator::Pointer duplicator = itk::ImageDuplicator::New(); duplicator->SetInputImage(itkImg); duplicator->Update(); parameters.m_SignalGen.m_FrequencyMap = duplicator->GetOutput(); parameters.m_Misc.m_ArtifactModelString += "_DISTORTED"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Distortions", BoolProperty::New(true)); } } parameters.m_SignalGen.m_EddyStrength = 0; parameters.m_Misc.m_CheckAddEddyCurrentsBox = m_Controls->m_AddEddy->isChecked(); if (m_Controls->m_AddEddy->isChecked()) { parameters.m_SignalGen.m_EddyStrength = m_Controls->m_EddyGradientStrength->value(); parameters.m_Misc.m_ArtifactModelString += "_EDDY"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Eddy-strength", DoubleProperty::New(parameters.m_SignalGen.m_EddyStrength)); } // Motion parameters.m_SignalGen.m_DoAddMotion = m_Controls->m_AddMotion->isChecked(); parameters.m_SignalGen.m_DoRandomizeMotion = m_Controls->m_RandomMotion->isChecked(); parameters.m_SignalGen.m_Translation[0] = m_Controls->m_MaxTranslationBoxX->value(); parameters.m_SignalGen.m_Translation[1] = m_Controls->m_MaxTranslationBoxY->value(); parameters.m_SignalGen.m_Translation[2] = m_Controls->m_MaxTranslationBoxZ->value(); parameters.m_SignalGen.m_Rotation[0] = m_Controls->m_MaxRotationBoxX->value(); parameters.m_SignalGen.m_Rotation[1] = m_Controls->m_MaxRotationBoxY->value(); parameters.m_SignalGen.m_Rotation[2] = m_Controls->m_MaxRotationBoxZ->value(); if ( m_Controls->m_AddMotion->isChecked() && m_Controls->m_FiberBundleComboBox->GetSelectedNode().IsNotNull() ) { parameters.m_Misc.m_ArtifactModelString += "_MOTION"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Random", BoolProperty::New(parameters.m_SignalGen.m_DoRandomizeMotion)); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Translation-x", DoubleProperty::New(parameters.m_SignalGen.m_Translation[0])); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Translation-y", DoubleProperty::New(parameters.m_SignalGen.m_Translation[1])); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Translation-z", DoubleProperty::New(parameters.m_SignalGen.m_Translation[2])); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Rotation-x", DoubleProperty::New(parameters.m_SignalGen.m_Rotation[0])); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Rotation-y", DoubleProperty::New(parameters.m_SignalGen.m_Rotation[1])); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Rotation-z", DoubleProperty::New(parameters.m_SignalGen.m_Rotation[2])); } // other imaging parameters + parameters.m_SignalGen.m_PartialFourier = m_Controls->m_PartialFourier->value(); parameters.m_SignalGen.m_ReversePhase = m_Controls->m_ReversePhaseBox->isChecked(); parameters.m_SignalGen.m_tLine = m_Controls->m_LineReadoutTimeBox->value(); parameters.m_SignalGen.m_tInhom = m_Controls->m_T2starBox->value(); parameters.m_SignalGen.m_tEcho = m_Controls->m_TEbox->value(); parameters.m_SignalGen.m_tRep = m_Controls->m_TRbox->value(); parameters.m_SignalGen.m_DoDisablePartialVolume = m_Controls->m_EnforcePureFiberVoxelsBox->isChecked(); parameters.m_SignalGen.m_AxonRadius = m_Controls->m_FiberRadius->value(); parameters.m_SignalGen.m_SignalScale = m_Controls->m_SignalScaleBox->value(); // adjust echo time if needed if ( parameters.m_SignalGen.m_tEcho < parameters.m_SignalGen.m_ImageRegion.GetSize(1)*parameters.m_SignalGen.m_tLine ) { this->m_Controls->m_TEbox->setValue( parameters.m_SignalGen.m_ImageRegion.GetSize(1)*parameters.m_SignalGen.m_tLine ); parameters.m_SignalGen.m_tEcho = m_Controls->m_TEbox->value(); QMessageBox::information( NULL, "Warning", "Echo time is too short! Time not sufficient to read slice. Automaticall adjusted to "+QString::number(parameters.m_SignalGen.m_tEcho)+" ms"); } // Noise parameters.m_Misc.m_CheckAddNoiseBox = m_Controls->m_AddNoise->isChecked(); + parameters.m_SignalGen.m_NoiseVariance = 0; if (m_Controls->m_AddNoise->isChecked()) { double noiseVariance = m_Controls->m_NoiseLevel->value(); + + // { + // switch (m_Controls->m_NoiseDistributionBox->currentIndex()) + // { + // case 0: + // { + // parameters.m_NoiseModel = new mitk::RicianNoiseModel(); + // parameters.m_Misc.m_ArtifactModelString += "_RICIAN-"; + // parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Noise-Distribution", StringProperty::New("Rician")); + // break; + // } + // case 1: + // { + // parameters.m_NoiseModel = new mitk::ChiSquareNoiseModel(); + // parameters.m_Misc.m_ArtifactModelString += "_CHISQUARED-"; + // parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Noise-Distribution", StringProperty::New("Chi-squared")); + // break; + // } + // default: + // { + // parameters.m_NoiseModel = new mitk::RicianNoiseModel(); + // parameters.m_Misc.m_ArtifactModelString += "_RICIAN-"; + // parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Noise-Distribution", StringProperty::New("Rician")); + // } + // } + // } + // parameters.m_NoiseModel->SetNoiseVariance(noiseVariance); + parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; + if (noiseVariance>0) { - switch (m_Controls->m_NoiseDistributionBox->currentIndex()) - { - case 0: - { - parameters.m_NoiseModel = new mitk::RicianNoiseModel(); - parameters.m_Misc.m_ArtifactModelString += "_RICIAN-"; - parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Noise-Distribution", StringProperty::New("Rician")); - break; - } - case 1: - { - parameters.m_NoiseModel = new mitk::ChiSquareNoiseModel(); - parameters.m_Misc.m_ArtifactModelString += "_CHISQUARED-"; - parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Noise-Distribution", StringProperty::New("Chi-squared")); - break; - } - default: - { - parameters.m_NoiseModel = new mitk::RicianNoiseModel(); - parameters.m_Misc.m_ArtifactModelString += "_RICIAN-"; - parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Noise-Distribution", StringProperty::New("Rician")); - } - } + parameters.m_SignalGen.m_NoiseVariance = m_Controls->m_NoiseLevel->value(); + parameters.m_Misc.m_ArtifactModelString += QString::number(noiseVariance).toStdString(); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Noise-Variance", DoubleProperty::New(noiseVariance)); } - parameters.m_NoiseModel->SetNoiseVariance(noiseVariance); - parameters.m_Misc.m_ArtifactModelString += QString::number(noiseVariance).toStdString(); - parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Noise-Variance", DoubleProperty::New(noiseVariance)); } // adjusting line readout time to the adapted image size needed for the DFT unsigned int y = parameters.m_SignalGen.m_ImageRegion.GetSize(1); y += y%2; if ( y>parameters.m_SignalGen.m_ImageRegion.GetSize(1) ) parameters.m_SignalGen.m_tLine *= (double)parameters.m_SignalGen.m_ImageRegion.GetSize(1)/y; // signal models { // compartment 1 switch (m_Controls->m_Compartment1Box->currentIndex()) { case 0: { mitk::StickModel* model = new mitk::StickModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(parameters.m_SignalGen.m_Bvalue); model->SetDiffusivity(m_Controls->m_StickWidget1->GetD()); model->SetT2(m_Controls->m_StickWidget1->GetT2()); model->SetT1(m_Controls->m_StickWidget1->GetT1()); model->m_CompartmentId = 1; parameters.m_FiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Stick"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Stick") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D", DoubleProperty::New(m_Controls->m_StickWidget1->GetD()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.T2", DoubleProperty::New(model->GetT2()) ); break; } case 1: { mitk::TensorModel* model = new mitk::TensorModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(parameters.m_SignalGen.m_Bvalue); model->SetDiffusivity1(m_Controls->m_ZeppelinWidget1->GetD1()); model->SetDiffusivity2(m_Controls->m_ZeppelinWidget1->GetD2()); model->SetDiffusivity3(m_Controls->m_ZeppelinWidget1->GetD2()); model->SetT2(m_Controls->m_ZeppelinWidget1->GetT2()); model->SetT1(m_Controls->m_ZeppelinWidget1->GetT1()); model->m_CompartmentId = 1; parameters.m_FiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Zeppelin"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Zeppelin") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D1", DoubleProperty::New(m_Controls->m_ZeppelinWidget1->GetD1()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D2", DoubleProperty::New(m_Controls->m_ZeppelinWidget1->GetD2()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.T2", DoubleProperty::New(model->GetT2()) ); break; } case 2: { mitk::TensorModel* model = new mitk::TensorModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(parameters.m_SignalGen.m_Bvalue); model->SetDiffusivity1(m_Controls->m_TensorWidget1->GetD1()); model->SetDiffusivity2(m_Controls->m_TensorWidget1->GetD2()); model->SetDiffusivity3(m_Controls->m_TensorWidget1->GetD3()); model->SetT2(m_Controls->m_TensorWidget1->GetT2()); model->SetT1(m_Controls->m_TensorWidget1->GetT1()); model->m_CompartmentId = 1; parameters.m_FiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Tensor"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Tensor") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D1", DoubleProperty::New(m_Controls->m_TensorWidget1->GetD1()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D2", DoubleProperty::New(m_Controls->m_TensorWidget1->GetD2()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D3", DoubleProperty::New(m_Controls->m_TensorWidget1->GetD3()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.T2", DoubleProperty::New(model->GetT2()) ); break; } case 3: { mitk::RawShModel* model = new mitk::RawShModel(); parameters.m_SignalGen.m_DoSimulateRelaxation = false; model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetMaxNumKernels(m_Controls->m_PrototypeWidget1->GetNumberOfSamples()); model->SetFaRange(m_Controls->m_PrototypeWidget1->GetMinFa(), m_Controls->m_PrototypeWidget1->GetMaxFa()); model->SetAdcRange(m_Controls->m_PrototypeWidget1->GetMinAdc(), m_Controls->m_PrototypeWidget1->GetMaxAdc()); model->m_CompartmentId = 1; parameters.m_FiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Prototype"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Prototype") ); break; } } if (m_Controls->m_Comp1VolumeFraction->GetSelectedNode().IsNotNull()) { mitk::DataNode::Pointer volumeNode = m_Controls->m_Comp1VolumeFraction->GetSelectedNode(); ItkDoubleImgType::Pointer comp1VolumeImage = ItkDoubleImgType::New(); mitk::Image* img = dynamic_cast(volumeNode->GetData()); CastToItkImage< ItkDoubleImgType >(img, comp1VolumeImage); parameters.m_FiberModelList.back()->SetVolumeFractionImage(comp1VolumeImage); } // compartment 2 switch (m_Controls->m_Compartment2Box->currentIndex()) { case 0: break; case 1: { mitk::StickModel* model = new mitk::StickModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(parameters.m_SignalGen.m_Bvalue); model->SetDiffusivity(m_Controls->m_StickWidget2->GetD()); model->SetT2(m_Controls->m_StickWidget2->GetT2()); model->SetT1(m_Controls->m_StickWidget2->GetT1()); model->m_CompartmentId = 2; parameters.m_FiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Stick"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Description", StringProperty::New("Inter-axonal compartment") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Model", StringProperty::New("Stick") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D", DoubleProperty::New(m_Controls->m_StickWidget2->GetD()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.T2", DoubleProperty::New(model->GetT2()) ); break; } case 2: { mitk::TensorModel* model = new mitk::TensorModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(parameters.m_SignalGen.m_Bvalue); model->SetDiffusivity1(m_Controls->m_ZeppelinWidget2->GetD1()); model->SetDiffusivity2(m_Controls->m_ZeppelinWidget2->GetD2()); model->SetDiffusivity3(m_Controls->m_ZeppelinWidget2->GetD2()); model->SetT2(m_Controls->m_ZeppelinWidget2->GetT2()); model->SetT1(m_Controls->m_ZeppelinWidget2->GetT1()); model->m_CompartmentId = 2; parameters.m_FiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Zeppelin"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Description", StringProperty::New("Inter-axonal compartment") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Model", StringProperty::New("Zeppelin") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D1", DoubleProperty::New(m_Controls->m_ZeppelinWidget2->GetD1()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D2", DoubleProperty::New(m_Controls->m_ZeppelinWidget2->GetD2()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.T2", DoubleProperty::New(model->GetT2()) ); break; } case 3: { mitk::TensorModel* model = new mitk::TensorModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(parameters.m_SignalGen.m_Bvalue); model->SetDiffusivity1(m_Controls->m_TensorWidget2->GetD1()); model->SetDiffusivity2(m_Controls->m_TensorWidget2->GetD2()); model->SetDiffusivity3(m_Controls->m_TensorWidget2->GetD3()); model->SetT2(m_Controls->m_TensorWidget2->GetT2()); model->SetT1(m_Controls->m_TensorWidget2->GetT1()); model->m_CompartmentId = 2; parameters.m_FiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Tensor"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Description", StringProperty::New("Inter-axonal compartment") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Model", StringProperty::New("Tensor") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D1", DoubleProperty::New(m_Controls->m_TensorWidget2->GetD1()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D2", DoubleProperty::New(m_Controls->m_TensorWidget2->GetD2()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D3", DoubleProperty::New(m_Controls->m_TensorWidget2->GetD3()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.T2", DoubleProperty::New(model->GetT2()) ); break; } } if (m_Controls->m_Comp2VolumeFraction->GetSelectedNode().IsNotNull() && parameters.m_FiberModelList.size()==2) { mitk::DataNode::Pointer volumeNode = m_Controls->m_Comp2VolumeFraction->GetSelectedNode(); ItkDoubleImgType::Pointer comp1VolumeImage = ItkDoubleImgType::New(); mitk::Image* img = dynamic_cast(volumeNode->GetData()); CastToItkImage< ItkDoubleImgType >(img, comp1VolumeImage); parameters.m_FiberModelList.back()->SetVolumeFractionImage(comp1VolumeImage); } // compartment 3 switch (m_Controls->m_Compartment3Box->currentIndex()) { case 0: { mitk::BallModel* model = new mitk::BallModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(parameters.m_SignalGen.m_Bvalue); model->SetDiffusivity(m_Controls->m_BallWidget1->GetD()); model->SetT2(m_Controls->m_BallWidget1->GetT2()); model->SetT1(m_Controls->m_BallWidget1->GetT1()); model->m_CompartmentId = 3; parameters.m_NonFiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Ball"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Ball") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.D", DoubleProperty::New(m_Controls->m_BallWidget1->GetD()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.T2", DoubleProperty::New(model->GetT2()) ); break; } case 1: { mitk::AstroStickModel* model = new mitk::AstroStickModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(parameters.m_SignalGen.m_Bvalue); model->SetDiffusivity(m_Controls->m_AstrosticksWidget1->GetD()); model->SetT2(m_Controls->m_AstrosticksWidget1->GetT2()); model->SetT1(m_Controls->m_AstrosticksWidget1->GetT1()); model->SetRandomizeSticks(m_Controls->m_AstrosticksWidget1->GetRandomizeSticks()); model->m_CompartmentId = 3; parameters.m_NonFiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Astrosticks"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Astrosticks") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.D", DoubleProperty::New(m_Controls->m_AstrosticksWidget1->GetD()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.T2", DoubleProperty::New(model->GetT2()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.RandomSticks", BoolProperty::New(m_Controls->m_AstrosticksWidget1->GetRandomizeSticks()) ); break; } case 2: { mitk::DotModel* model = new mitk::DotModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetT2(m_Controls->m_DotWidget1->GetT2()); model->SetT1(m_Controls->m_DotWidget1->GetT1()); model->m_CompartmentId = 3; parameters.m_NonFiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Dot"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Dot") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.T2", DoubleProperty::New(model->GetT2()) ); break; } case 3: { mitk::RawShModel* model = new mitk::RawShModel(); parameters.m_SignalGen.m_DoSimulateRelaxation = false; model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetMaxNumKernels(m_Controls->m_PrototypeWidget3->GetNumberOfSamples()); model->SetFaRange(m_Controls->m_PrototypeWidget3->GetMinFa(), m_Controls->m_PrototypeWidget3->GetMaxFa()); model->SetAdcRange(m_Controls->m_PrototypeWidget3->GetMinAdc(), m_Controls->m_PrototypeWidget3->GetMaxAdc()); model->m_CompartmentId = 3; parameters.m_NonFiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Prototype"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Prototype") ); break; } } if (m_Controls->m_Comp3VolumeFraction->GetSelectedNode().IsNotNull()) { mitk::DataNode::Pointer volumeNode = m_Controls->m_Comp3VolumeFraction->GetSelectedNode(); ItkDoubleImgType::Pointer comp1VolumeImage = ItkDoubleImgType::New(); mitk::Image* img = dynamic_cast(volumeNode->GetData()); CastToItkImage< ItkDoubleImgType >(img, comp1VolumeImage); parameters.m_NonFiberModelList.back()->SetVolumeFractionImage(comp1VolumeImage); } switch (m_Controls->m_Compartment4Box->currentIndex()) { case 0: break; case 1: { mitk::BallModel* model = new mitk::BallModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(parameters.m_SignalGen.m_Bvalue); model->SetDiffusivity(m_Controls->m_BallWidget2->GetD()); model->SetT2(m_Controls->m_BallWidget2->GetT2()); model->SetT1(m_Controls->m_BallWidget2->GetT1()); model->m_CompartmentId = 4; parameters.m_NonFiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Ball"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Ball") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.D", DoubleProperty::New(m_Controls->m_BallWidget2->GetD()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.T2", DoubleProperty::New(model->GetT2()) ); break; } case 2: { mitk::AstroStickModel* model = new mitk::AstroStickModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(parameters.m_SignalGen.m_Bvalue); model->SetDiffusivity(m_Controls->m_AstrosticksWidget2->GetD()); model->SetT2(m_Controls->m_AstrosticksWidget2->GetT2()); model->SetT1(m_Controls->m_AstrosticksWidget2->GetT1()); model->SetRandomizeSticks(m_Controls->m_AstrosticksWidget2->GetRandomizeSticks()); model->m_CompartmentId = 4; parameters.m_NonFiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Astrosticks"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Astrosticks") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.D", DoubleProperty::New(m_Controls->m_AstrosticksWidget2->GetD()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.T2", DoubleProperty::New(model->GetT2()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.RandomSticks", BoolProperty::New(m_Controls->m_AstrosticksWidget2->GetRandomizeSticks()) ); break; } case 3: { mitk::DotModel* model = new mitk::DotModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetT2(m_Controls->m_DotWidget2->GetT2()); model->SetT1(m_Controls->m_DotWidget2->GetT1()); model->m_CompartmentId = 4; parameters.m_NonFiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Dot"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Dot") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.T2", DoubleProperty::New(model->GetT2()) ); break; } case 4: { mitk::RawShModel* model = new mitk::RawShModel(); parameters.m_SignalGen.m_DoSimulateRelaxation = false; model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetMaxNumKernels(m_Controls->m_PrototypeWidget4->GetNumberOfSamples()); model->SetFaRange(m_Controls->m_PrototypeWidget4->GetMinFa(), m_Controls->m_PrototypeWidget4->GetMaxFa()); model->SetAdcRange(m_Controls->m_PrototypeWidget4->GetMinAdc(), m_Controls->m_PrototypeWidget4->GetMaxAdc()); model->m_CompartmentId = 4; parameters.m_NonFiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Prototype"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Prototype") ); break; } } if (m_Controls->m_Comp4VolumeFraction->GetSelectedNode().IsNotNull() && parameters.m_NonFiberModelList.size()==2) { mitk::DataNode::Pointer volumeNode = m_Controls->m_Comp4VolumeFraction->GetSelectedNode(); ItkDoubleImgType::Pointer compVolumeImage = ItkDoubleImgType::New(); mitk::Image* img = dynamic_cast(volumeNode->GetData()); CastToItkImage< ItkDoubleImgType >(img, compVolumeImage); parameters.m_NonFiberModelList.back()->SetVolumeFractionImage(compVolumeImage); } } // check if comp 3 or 4 volume fraction image is set if (parameters.m_NonFiberModelList.size()==2 && (parameters.m_NonFiberModelList[0]->GetVolumeFractionImage()==nullptr || parameters.m_NonFiberModelList[1]->GetVolumeFractionImage()==nullptr)) { m_Controls->m_Compartment4Box->setCurrentIndex(0); parameters.m_NonFiberModelList.pop_back(); QMessageBox::information(NULL, "Compartment 4 disabled", "More than one non-fiber compartment selected but no volume fraction maps set!"); } -// RELICT -// parameters.m_SignalGen.m_FiberSeparationThreshold = m_Controls->m_SeparationAngleBox->value(); -// switch (m_Controls->m_DiffusionDirectionBox->currentIndex()) -// { -// case 0: -// parameters.m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::FIBER_TANGENT_DIRECTIONS; -// break; -// case 1: -// parameters.m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::MAIN_FIBER_DIRECTIONS; -// break; -// case 2: -// parameters.m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::RANDOM_DIRECTIONS; -// parameters.m_SignalGen.m_DoAddMotion = false; -// parameters.m_SignalGen.m_DoAddGibbsRinging = false; -// parameters.m_SignalGen.m_KspaceLineOffset = 0.0; -// parameters.m_SignalGen.m_FrequencyMap = NULL; -// parameters.m_SignalGen.m_CroppingFactor = 1.0; -// parameters.m_SignalGen.m_EddyStrength = 0; -// break; -// default: -// parameters.m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::FIBER_TANGENT_DIRECTIONS; -// } + // RELICT + // parameters.m_SignalGen.m_FiberSeparationThreshold = m_Controls->m_SeparationAngleBox->value(); + // switch (m_Controls->m_DiffusionDirectionBox->currentIndex()) + // { + // case 0: + // parameters.m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::FIBER_TANGENT_DIRECTIONS; + // break; + // case 1: + // parameters.m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::MAIN_FIBER_DIRECTIONS; + // break; + // case 2: + // parameters.m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::RANDOM_DIRECTIONS; + // parameters.m_SignalGen.m_DoAddMotion = false; + // parameters.m_SignalGen.m_DoAddGibbsRinging = false; + // parameters.m_SignalGen.m_KspaceLineOffset = 0.0; + // parameters.m_SignalGen.m_FrequencyMap = NULL; + // parameters.m_SignalGen.m_CroppingFactor = 1.0; + // parameters.m_SignalGen.m_EddyStrength = 0; + // break; + // default: + // parameters.m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::FIBER_TANGENT_DIRECTIONS; + // } parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.SignalScale", IntProperty::New(parameters.m_SignalGen.m_SignalScale)); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.FiberRadius", IntProperty::New(parameters.m_SignalGen.m_AxonRadius)); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Tinhom", DoubleProperty::New(parameters.m_SignalGen.m_tInhom)); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Tline", DoubleProperty::New(parameters.m_SignalGen.m_tLine)); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.TE", DoubleProperty::New(parameters.m_SignalGen.m_tEcho)); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.b-value", DoubleProperty::New(parameters.m_SignalGen.m_Bvalue)); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.NoPartialVolume", BoolProperty::New(parameters.m_SignalGen.m_DoDisablePartialVolume)); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Relaxation", BoolProperty::New(parameters.m_SignalGen.m_DoSimulateRelaxation)); parameters.m_Misc.m_ResultNode->AddProperty("binary", BoolProperty::New(false)); parameters.m_Misc.m_CheckRealTimeFibersBox = m_Controls->m_RealTimeFibers->isChecked(); parameters.m_Misc.m_CheckAdvancedFiberOptionsBox = m_Controls->m_AdvancedOptionsBox->isChecked(); parameters.m_Misc.m_CheckIncludeFiducialsBox = m_Controls->m_IncludeFiducials->isChecked(); parameters.m_Misc.m_CheckConstantRadiusBox = m_Controls->m_ConstantRadiusBox->isChecked(); switch(m_Controls->m_DistributionBox->currentIndex()) { case 0: parameters.m_FiberGen.m_Distribution = FiberGenerationParameters::DISTRIBUTE_UNIFORM; break; case 1: parameters.m_FiberGen.m_Distribution = FiberGenerationParameters::DISTRIBUTE_GAUSSIAN; break; default: parameters.m_FiberGen.m_Distribution = FiberGenerationParameters::DISTRIBUTE_UNIFORM; } parameters.m_FiberGen.m_Variance = m_Controls->m_VarianceBox->value(); parameters.m_FiberGen.m_Density = m_Controls->m_FiberDensityBox->value(); parameters.m_FiberGen.m_Sampling = m_Controls->m_FiberSamplingBox->value(); parameters.m_FiberGen.m_Tension = m_Controls->m_TensionBox->value(); parameters.m_FiberGen.m_Continuity = m_Controls->m_ContinuityBox->value(); parameters.m_FiberGen.m_Bias = m_Controls->m_BiasBox->value(); parameters.m_FiberGen.m_Rotation[0] = m_Controls->m_XrotBox->value(); parameters.m_FiberGen.m_Rotation[1] = m_Controls->m_YrotBox->value(); parameters.m_FiberGen.m_Rotation[2] = m_Controls->m_ZrotBox->value(); parameters.m_FiberGen.m_Translation[0] = m_Controls->m_XtransBox->value(); parameters.m_FiberGen.m_Translation[1] = m_Controls->m_YtransBox->value(); parameters.m_FiberGen.m_Translation[2] = m_Controls->m_ZtransBox->value(); parameters.m_FiberGen.m_Scale[0] = m_Controls->m_XscaleBox->value(); parameters.m_FiberGen.m_Scale[1] = m_Controls->m_YscaleBox->value(); parameters.m_FiberGen.m_Scale[2] = m_Controls->m_ZscaleBox->value(); return parameters; } void QmitkFiberfoxView::SaveParameters() { FiberfoxParameters<> ffParamaters = UpdateImageParameters(); QString filename = QFileDialog::getSaveFileName( 0, tr("Save Parameters"), m_ParameterFile, tr("Fiberfox Parameters (*.ffp)") ); bool ok = true; bool first = true; bool dosampling = false; mitk::Image::Pointer diffImg = NULL; itk::Image< itk::DiffusionTensor3D< double >, 3 >::Pointer tensorImage = NULL; const int shOrder = 2; typedef itk::AnalyticalDiffusionQballReconstructionImageFilter QballFilterType; QballFilterType::CoefficientImageType::Pointer itkFeatureImage = NULL; ItkDoubleImgType::Pointer adcImage = NULL; for (unsigned int i=0; i* model = NULL; if (i* >(ffParamaters.m_FiberModelList.at(i)); else model = dynamic_cast< mitk::RawShModel<>* >(ffParamaters.m_NonFiberModelList.at(i-ffParamaters.m_FiberModelList.size())); if (model!=0 && model->GetNumberOfKernels()<=0) { if (first==true) { if (QMessageBox::question(NULL, "Prototype signal sampling", "Do you want to sample prototype signals from the selected diffusion-weighted imag and save them?",QMessageBox::Yes,QMessageBox::No)==QMessageBox::Yes) dosampling = true; first = false; if (dosampling && (m_Controls->m_TemplateComboBox->GetSelectedNode().IsNull() || !mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_Controls->m_TemplateComboBox->GetSelectedNode()->GetData())))) { QMessageBox::information(NULL, "Parameter file not saved", "No diffusion-weighted image selected to sample signal from."); return; } else if (dosampling) { diffImg = dynamic_cast(m_Controls->m_TemplateComboBox->GetSelectedNode()->GetData()); typedef itk::DiffusionTensor3DReconstructionImageFilter< short, short, double > TensorReconstructionImageFilterType; TensorReconstructionImageFilterType::Pointer filter = TensorReconstructionImageFilterType::New(); ItkDwiType::Pointer itkVectorImagePointer = ItkDwiType::New(); mitk::CastToItkImage(diffImg, itkVectorImagePointer); filter->SetGradientImage( static_cast( diffImg->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer(), itkVectorImagePointer ); filter->SetBValue( static_cast(diffImg->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue() ); filter->Update(); tensorImage = filter->GetOutput(); const int NumCoeffs = (shOrder*shOrder + shOrder + 2)/2 + shOrder; QballFilterType::Pointer qballfilter = QballFilterType::New(); qballfilter->SetGradientImage( static_cast( diffImg->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer(), itkVectorImagePointer ); qballfilter->SetBValue( static_cast(diffImg->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue() ); qballfilter->SetLambda(0.006); qballfilter->SetNormalizationMethod(QballFilterType::QBAR_RAW_SIGNAL); qballfilter->Update(); itkFeatureImage = qballfilter->GetCoefficientImage(); itk::AdcImageFilter< short, double >::Pointer adcFilter = itk::AdcImageFilter< short, double >::New(); adcFilter->SetInput( itkVectorImagePointer ); adcFilter->SetGradientDirections( static_cast( diffImg->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer() ); adcFilter->SetB_value( static_cast(diffImg->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue() ); adcFilter->Update(); adcImage = adcFilter->GetOutput(); } } if (dosampling && diffImg.IsNotNull()) { ok = model->SampleKernels(diffImg, ffParamaters.m_SignalGen.m_MaskImage, tensorImage, itkFeatureImage, adcImage); if (!ok) { QMessageBox::information( NULL, "Parameter file not saved", "No valid prototype signals could be sampled."); return; } } } } ffParamaters.SaveParameters(filename.toStdString()); m_ParameterFile = filename; } void QmitkFiberfoxView::LoadParameters() { QString filename = QFileDialog::getOpenFileName(0, tr("Load Parameters"), QString(itksys::SystemTools::GetFilenamePath(m_ParameterFile.toStdString()).c_str()), tr("Fiberfox Parameters (*.ffp)") ); if(filename.isEmpty() || filename.isNull()) return; m_ParameterFile = filename; FiberfoxParameters<> parameters; parameters.LoadParameters(filename.toStdString()); m_Controls->m_RealTimeFibers->setChecked(parameters.m_Misc.m_CheckRealTimeFibersBox); m_Controls->m_AdvancedOptionsBox->setChecked(parameters.m_Misc.m_CheckAdvancedFiberOptionsBox); m_Controls->m_IncludeFiducials->setChecked(parameters.m_Misc.m_CheckIncludeFiducialsBox); m_Controls->m_ConstantRadiusBox->setChecked(parameters.m_Misc.m_CheckConstantRadiusBox); m_Controls->m_DistributionBox->setCurrentIndex(parameters.m_FiberGen.m_Distribution); m_Controls->m_VarianceBox->setValue(parameters.m_FiberGen.m_Variance); m_Controls->m_FiberDensityBox->setValue(parameters.m_FiberGen.m_Density); m_Controls->m_FiberSamplingBox->setValue(parameters.m_FiberGen.m_Sampling); m_Controls->m_TensionBox->setValue(parameters.m_FiberGen.m_Tension); m_Controls->m_ContinuityBox->setValue(parameters.m_FiberGen.m_Continuity); m_Controls->m_BiasBox->setValue(parameters.m_FiberGen.m_Bias); m_Controls->m_XrotBox->setValue(parameters.m_FiberGen.m_Rotation[0]); m_Controls->m_YrotBox->setValue(parameters.m_FiberGen.m_Rotation[1]); m_Controls->m_ZrotBox->setValue(parameters.m_FiberGen.m_Rotation[2]); m_Controls->m_XtransBox->setValue(parameters.m_FiberGen.m_Translation[0]); m_Controls->m_YtransBox->setValue(parameters.m_FiberGen.m_Translation[1]); m_Controls->m_ZtransBox->setValue(parameters.m_FiberGen.m_Translation[2]); m_Controls->m_XscaleBox->setValue(parameters.m_FiberGen.m_Scale[0]); m_Controls->m_YscaleBox->setValue(parameters.m_FiberGen.m_Scale[1]); m_Controls->m_ZscaleBox->setValue(parameters.m_FiberGen.m_Scale[2]); // image generation parameters m_Controls->m_SizeX->setValue(parameters.m_SignalGen.m_ImageRegion.GetSize(0)); m_Controls->m_SizeY->setValue(parameters.m_SignalGen.m_ImageRegion.GetSize(1)); m_Controls->m_SizeZ->setValue(parameters.m_SignalGen.m_ImageRegion.GetSize(2)); m_Controls->m_SpacingX->setValue(parameters.m_SignalGen.m_ImageSpacing[0]); m_Controls->m_SpacingY->setValue(parameters.m_SignalGen.m_ImageSpacing[1]); m_Controls->m_SpacingZ->setValue(parameters.m_SignalGen.m_ImageSpacing[2]); m_Controls->m_NumGradientsBox->setValue(parameters.m_SignalGen.GetNumWeightedVolumes()); m_Controls->m_BvalueBox->setValue(parameters.m_SignalGen.m_Bvalue); m_Controls->m_SignalScaleBox->setValue(parameters.m_SignalGen.m_SignalScale); m_Controls->m_TEbox->setValue(parameters.m_SignalGen.m_tEcho); m_Controls->m_LineReadoutTimeBox->setValue(parameters.m_SignalGen.m_tLine); m_Controls->m_T2starBox->setValue(parameters.m_SignalGen.m_tInhom); m_Controls->m_FiberRadius->setValue(parameters.m_SignalGen.m_AxonRadius); m_Controls->m_RelaxationBox->setChecked(parameters.m_SignalGen.m_DoSimulateRelaxation); m_Controls->m_EnforcePureFiberVoxelsBox->setChecked(parameters.m_SignalGen.m_DoDisablePartialVolume); if (parameters.m_NoiseModel!=NULL) { m_Controls->m_AddNoise->setChecked(parameters.m_Misc.m_CheckAddNoiseBox); if (dynamic_cast*>(parameters.m_NoiseModel)) m_Controls->m_NoiseDistributionBox->setCurrentIndex(0); else if (dynamic_cast*>(parameters.m_NoiseModel)) m_Controls->m_NoiseDistributionBox->setCurrentIndex(1); m_Controls->m_NoiseLevel->setValue(parameters.m_NoiseModel->GetNoiseVariance()); } else m_Controls->m_AddNoise->setChecked(false); m_Controls->m_VolumeFractionsBox->setChecked(parameters.m_Misc.m_CheckOutputVolumeFractionsBox); m_Controls->m_AdvancedOptionsBox_2->setChecked(parameters.m_Misc.m_CheckAdvancedSignalOptionsBox); m_Controls->m_AddGhosts->setChecked(parameters.m_Misc.m_CheckAddGhostsBox); m_Controls->m_AddAliasing->setChecked(parameters.m_Misc.m_CheckAddAliasingBox); m_Controls->m_AddDistortions->setChecked(parameters.m_Misc.m_CheckAddDistortionsBox); m_Controls->m_AddSpikes->setChecked(parameters.m_Misc.m_CheckAddSpikesBox); m_Controls->m_AddEddy->setChecked(parameters.m_Misc.m_CheckAddEddyCurrentsBox); m_Controls->m_kOffsetBox->setValue(parameters.m_SignalGen.m_KspaceLineOffset); m_Controls->m_WrapBox->setValue(100*(1-parameters.m_SignalGen.m_CroppingFactor)); m_Controls->m_SpikeNumBox->setValue(parameters.m_SignalGen.m_Spikes); m_Controls->m_SpikeScaleBox->setValue(parameters.m_SignalGen.m_SpikeAmplitude); m_Controls->m_EddyGradientStrength->setValue(parameters.m_SignalGen.m_EddyStrength); m_Controls->m_AddGibbsRinging->setChecked(parameters.m_SignalGen.m_DoAddGibbsRinging); m_Controls->m_AddMotion->setChecked(parameters.m_SignalGen.m_DoAddMotion); m_Controls->m_RandomMotion->setChecked(parameters.m_SignalGen.m_DoRandomizeMotion); m_Controls->m_MaxTranslationBoxX->setValue(parameters.m_SignalGen.m_Translation[0]); m_Controls->m_MaxTranslationBoxY->setValue(parameters.m_SignalGen.m_Translation[1]); m_Controls->m_MaxTranslationBoxZ->setValue(parameters.m_SignalGen.m_Translation[2]); m_Controls->m_MaxRotationBoxX->setValue(parameters.m_SignalGen.m_Rotation[0]); m_Controls->m_MaxRotationBoxY->setValue(parameters.m_SignalGen.m_Rotation[1]); m_Controls->m_MaxRotationBoxZ->setValue(parameters.m_SignalGen.m_Rotation[2]); - m_Controls->m_DiffusionDirectionBox->setCurrentIndex(parameters.m_SignalGen.m_DiffusionDirectionMode); - m_Controls->m_SeparationAngleBox->setValue(parameters.m_SignalGen.m_FiberSeparationThreshold); m_Controls->m_Compartment1Box->setCurrentIndex(0); m_Controls->m_Compartment2Box->setCurrentIndex(0); m_Controls->m_Compartment3Box->setCurrentIndex(0); m_Controls->m_Compartment4Box->setCurrentIndex(0); for (unsigned int i=0; i* signalModel = NULL; if (im_CompartmentId) { case 1: { if (dynamic_cast*>(signalModel)) { mitk::StickModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_StickWidget1->SetT2(model->GetT2()); m_Controls->m_StickWidget1->SetD(model->GetDiffusivity()); m_Controls->m_Compartment1Box->setCurrentIndex(0); break; } else if (dynamic_cast*>(signalModel)) { mitk::TensorModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_TensorWidget1->SetT2(model->GetT2()); m_Controls->m_TensorWidget1->SetD1(model->GetDiffusivity1()); m_Controls->m_TensorWidget1->SetD2(model->GetDiffusivity2()); m_Controls->m_TensorWidget1->SetD3(model->GetDiffusivity3()); m_Controls->m_Compartment1Box->setCurrentIndex(2); break; } else if (dynamic_cast*>(signalModel)) { mitk::RawShModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_PrototypeWidget1->SetNumberOfSamples(model->GetMaxNumKernels()); m_Controls->m_PrototypeWidget1->SetMinFa(model->GetFaRange().first); m_Controls->m_PrototypeWidget1->SetMaxFa(model->GetFaRange().second); m_Controls->m_PrototypeWidget1->SetMinAdc(model->GetAdcRange().first); m_Controls->m_PrototypeWidget1->SetMaxAdc(model->GetAdcRange().second); m_Controls->m_Compartment1Box->setCurrentIndex(3); break; } break; } case 2: { if (dynamic_cast*>(signalModel)) { mitk::StickModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_StickWidget2->SetT2(model->GetT2()); m_Controls->m_StickWidget2->SetD(model->GetDiffusivity()); m_Controls->m_Compartment2Box->setCurrentIndex(1); break; } else if (dynamic_cast*>(signalModel)) { mitk::TensorModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_TensorWidget2->SetT2(model->GetT2()); m_Controls->m_TensorWidget2->SetD1(model->GetDiffusivity1()); m_Controls->m_TensorWidget2->SetD2(model->GetDiffusivity2()); m_Controls->m_TensorWidget2->SetD3(model->GetDiffusivity3()); m_Controls->m_Compartment2Box->setCurrentIndex(3); break; } break; } case 3: { if (dynamic_cast*>(signalModel)) { mitk::BallModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_BallWidget1->SetT2(model->GetT2()); m_Controls->m_BallWidget1->SetD(model->GetDiffusivity()); m_Controls->m_Compartment3Box->setCurrentIndex(0); break; } else if (dynamic_cast*>(signalModel)) { mitk::AstroStickModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_AstrosticksWidget1->SetT2(model->GetT2()); m_Controls->m_AstrosticksWidget1->SetD(model->GetDiffusivity()); m_Controls->m_AstrosticksWidget1->SetRandomizeSticks(model->GetRandomizeSticks()); m_Controls->m_Compartment3Box->setCurrentIndex(1); break; } else if (dynamic_cast*>(signalModel)) { mitk::DotModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_DotWidget1->SetT2(model->GetT2()); m_Controls->m_Compartment3Box->setCurrentIndex(2); break; } else if (dynamic_cast*>(signalModel)) { mitk::RawShModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_PrototypeWidget3->SetNumberOfSamples(model->GetMaxNumKernels()); m_Controls->m_PrototypeWidget3->SetMinFa(model->GetFaRange().first); m_Controls->m_PrototypeWidget3->SetMaxFa(model->GetFaRange().second); m_Controls->m_PrototypeWidget3->SetMinAdc(model->GetAdcRange().first); m_Controls->m_PrototypeWidget3->SetMaxAdc(model->GetAdcRange().second); m_Controls->m_Compartment3Box->setCurrentIndex(3); break; } break; } case 4: { if (dynamic_cast*>(signalModel)) { mitk::BallModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_BallWidget2->SetT2(model->GetT2()); m_Controls->m_BallWidget2->SetD(model->GetDiffusivity()); m_Controls->m_Compartment4Box->setCurrentIndex(1); break; } else if (dynamic_cast*>(signalModel)) { mitk::AstroStickModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_AstrosticksWidget2->SetT2(model->GetT2()); m_Controls->m_AstrosticksWidget2->SetD(model->GetDiffusivity()); m_Controls->m_AstrosticksWidget2->SetRandomizeSticks(model->GetRandomizeSticks()); m_Controls->m_Compartment4Box->setCurrentIndex(2); break; } else if (dynamic_cast*>(signalModel)) { mitk::DotModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_DotWidget2->SetT2(model->GetT2()); m_Controls->m_Compartment4Box->setCurrentIndex(3); break; } else if (dynamic_cast*>(signalModel)) { mitk::RawShModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_PrototypeWidget4->SetNumberOfSamples(model->GetMaxNumKernels()); m_Controls->m_PrototypeWidget4->SetMinFa(model->GetFaRange().first); m_Controls->m_PrototypeWidget4->SetMaxFa(model->GetFaRange().second); m_Controls->m_PrototypeWidget4->SetMinAdc(model->GetAdcRange().first); m_Controls->m_PrototypeWidget4->SetMaxAdc(model->GetAdcRange().second); m_Controls->m_Compartment4Box->setCurrentIndex(4); break; } break; } } } } void QmitkFiberfoxView::ShowAdvancedOptions(int state) { if (state) { m_Controls->m_AdvancedFiberOptionsFrame->setVisible(true); m_Controls->m_AdvancedSignalOptionsFrame->setVisible(true); m_Controls->m_AdvancedOptionsBox->setChecked(true); m_Controls->m_AdvancedOptionsBox_2->setChecked(true); } else { m_Controls->m_AdvancedFiberOptionsFrame->setVisible(false); m_Controls->m_AdvancedSignalOptionsFrame->setVisible(false); m_Controls->m_AdvancedOptionsBox->setChecked(false); m_Controls->m_AdvancedOptionsBox_2->setChecked(false); } } void QmitkFiberfoxView::Comp1ModelFrameVisibility(int index) { m_Controls->m_StickWidget1->setVisible(false); m_Controls->m_ZeppelinWidget1->setVisible(false); m_Controls->m_TensorWidget1->setVisible(false); m_Controls->m_PrototypeWidget1->setVisible(false); switch (index) { case 0: m_Controls->m_StickWidget1->setVisible(true); break; case 1: m_Controls->m_ZeppelinWidget1->setVisible(true); break; case 2: m_Controls->m_TensorWidget1->setVisible(true); break; case 3: m_Controls->m_PrototypeWidget1->setVisible(true); break; } } void QmitkFiberfoxView::Comp2ModelFrameVisibility(int index) { m_Controls->m_StickWidget2->setVisible(false); m_Controls->m_ZeppelinWidget2->setVisible(false); m_Controls->m_TensorWidget2->setVisible(false); m_Controls->m_Comp2FractionFrame->setVisible(false); switch (index) { case 0: break; case 1: m_Controls->m_StickWidget2->setVisible(true); m_Controls->m_Comp2FractionFrame->setVisible(true); break; case 2: m_Controls->m_ZeppelinWidget2->setVisible(true); m_Controls->m_Comp2FractionFrame->setVisible(true); break; case 3: m_Controls->m_TensorWidget2->setVisible(true); m_Controls->m_Comp2FractionFrame->setVisible(true); break; } } void QmitkFiberfoxView::Comp3ModelFrameVisibility(int index) { m_Controls->m_BallWidget1->setVisible(false); m_Controls->m_AstrosticksWidget1->setVisible(false); m_Controls->m_DotWidget1->setVisible(false); m_Controls->m_PrototypeWidget3->setVisible(false); switch (index) { case 0: m_Controls->m_BallWidget1->setVisible(true); break; case 1: m_Controls->m_AstrosticksWidget1->setVisible(true); break; case 2: m_Controls->m_DotWidget1->setVisible(true); break; case 3: m_Controls->m_PrototypeWidget3->setVisible(true); break; } } void QmitkFiberfoxView::Comp4ModelFrameVisibility(int index) { m_Controls->m_BallWidget2->setVisible(false); m_Controls->m_AstrosticksWidget2->setVisible(false); m_Controls->m_DotWidget2->setVisible(false); m_Controls->m_PrototypeWidget4->setVisible(false); m_Controls->m_Comp4FractionFrame->setVisible(false); switch (index) { case 0: break; case 1: m_Controls->m_BallWidget2->setVisible(true); m_Controls->m_Comp4FractionFrame->setVisible(true); break; case 2: m_Controls->m_AstrosticksWidget2->setVisible(true); m_Controls->m_Comp4FractionFrame->setVisible(true); break; case 3: m_Controls->m_DotWidget2->setVisible(true); m_Controls->m_Comp4FractionFrame->setVisible(true); break; case 4: m_Controls->m_PrototypeWidget4->setVisible(true); m_Controls->m_Comp4FractionFrame->setVisible(true); break; } } void QmitkFiberfoxView::OnConstantRadius(int value) { if (value>0 && m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnAddMotion(int value) { if (value>0) m_Controls->m_MotionArtifactFrame->setVisible(true); else m_Controls->m_MotionArtifactFrame->setVisible(false); } void QmitkFiberfoxView::OnAddAliasing(int value) { if (value>0) m_Controls->m_AliasingFrame->setVisible(true); else m_Controls->m_AliasingFrame->setVisible(false); } void QmitkFiberfoxView::OnAddSpikes(int value) { if (value>0) m_Controls->m_SpikeFrame->setVisible(true); else m_Controls->m_SpikeFrame->setVisible(false); } void QmitkFiberfoxView::OnAddEddy(int value) { if (value>0) m_Controls->m_EddyFrame->setVisible(true); else m_Controls->m_EddyFrame->setVisible(false); } void QmitkFiberfoxView::OnAddDistortions(int value) { if (value>0) m_Controls->m_DistortionsFrame->setVisible(true); else m_Controls->m_DistortionsFrame->setVisible(false); } void QmitkFiberfoxView::OnAddGhosts(int value) { if (value>0) m_Controls->m_GhostFrame->setVisible(true); else m_Controls->m_GhostFrame->setVisible(false); } void QmitkFiberfoxView::OnAddNoise(int value) { if (value>0) m_Controls->m_NoiseFrame->setVisible(true); else m_Controls->m_NoiseFrame->setVisible(false); } void QmitkFiberfoxView::OnDistributionChanged(int value) { if (value==1) m_Controls->m_VarianceBox->setVisible(true); else m_Controls->m_VarianceBox->setVisible(false); if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnVarianceChanged(double) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnFiberDensityChanged(int) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnFiberSamplingChanged(double) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnTensionChanged(double) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnContinuityChanged(double) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnBiasChanged(double) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::AlignOnGrid() { for (unsigned int i=0; i(m_SelectedFiducials.at(i)->GetData()); mitk::Point3D wc0 = pe->GetWorldControlPoint(0); mitk::DataStorage::SetOfObjects::ConstPointer parentFibs = GetDataStorage()->GetSources(m_SelectedFiducials.at(i)); for( mitk::DataStorage::SetOfObjects::const_iterator it = parentFibs->begin(); it != parentFibs->end(); ++it ) { mitk::DataNode::Pointer pFibNode = *it; if ( pFibNode.IsNotNull() && dynamic_cast(pFibNode->GetData()) ) { mitk::DataStorage::SetOfObjects::ConstPointer parentImgs = GetDataStorage()->GetSources(pFibNode); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = parentImgs->begin(); it2 != parentImgs->end(); ++it2 ) { mitk::DataNode::Pointer pImgNode = *it2; if ( pImgNode.IsNotNull() && dynamic_cast(pImgNode->GetData()) ) { mitk::Image::Pointer img = dynamic_cast(pImgNode->GetData()); mitk::BaseGeometry::Pointer geom = img->GetGeometry(); itk::Index<3> idx; geom->WorldToIndex(wc0, idx); mitk::Point3D cIdx; cIdx[0]=idx[0]; cIdx[1]=idx[1]; cIdx[2]=idx[2]; mitk::Point3D world; geom->IndexToWorld(cIdx,world); mitk::Vector3D trans = world - wc0; pe->GetGeometry()->Translate(trans); break; } } break; } } } for(unsigned int i=0; iGetSources(fibNode); for( mitk::DataStorage::SetOfObjects::const_iterator it = sources->begin(); it != sources->end(); ++it ) { mitk::DataNode::Pointer imgNode = *it; if ( imgNode.IsNotNull() && dynamic_cast(imgNode->GetData()) ) { mitk::DataStorage::SetOfObjects::ConstPointer derivations = GetDataStorage()->GetDerivations(fibNode); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = derivations->begin(); it2 != derivations->end(); ++it2 ) { mitk::DataNode::Pointer fiducialNode = *it2; if ( fiducialNode.IsNotNull() && dynamic_cast(fiducialNode->GetData()) ) { mitk::PlanarEllipse::Pointer pe = dynamic_cast(fiducialNode->GetData()); mitk::Point3D wc0 = pe->GetWorldControlPoint(0); mitk::Image::Pointer img = dynamic_cast(imgNode->GetData()); mitk::BaseGeometry::Pointer geom = img->GetGeometry(); itk::Index<3> idx; geom->WorldToIndex(wc0, idx); mitk::Point3D cIdx; cIdx[0]=idx[0]; cIdx[1]=idx[1]; cIdx[2]=idx[2]; mitk::Point3D world; geom->IndexToWorld(cIdx,world); mitk::Vector3D trans = world - wc0; pe->GetGeometry()->Translate(trans); } } break; } } } for(unsigned int i=0; i(m_SelectedImages.at(i)->GetData()); mitk::DataStorage::SetOfObjects::ConstPointer derivations = GetDataStorage()->GetDerivations(m_SelectedImages.at(i)); for( mitk::DataStorage::SetOfObjects::const_iterator it = derivations->begin(); it != derivations->end(); ++it ) { mitk::DataNode::Pointer fibNode = *it; if ( fibNode.IsNotNull() && dynamic_cast(fibNode->GetData()) ) { mitk::DataStorage::SetOfObjects::ConstPointer derivations2 = GetDataStorage()->GetDerivations(fibNode); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = derivations2->begin(); it2 != derivations2->end(); ++it2 ) { mitk::DataNode::Pointer fiducialNode = *it2; if ( fiducialNode.IsNotNull() && dynamic_cast(fiducialNode->GetData()) ) { mitk::PlanarEllipse::Pointer pe = dynamic_cast(fiducialNode->GetData()); mitk::Point3D wc0 = pe->GetWorldControlPoint(0); mitk::BaseGeometry::Pointer geom = img->GetGeometry(); itk::Index<3> idx; geom->WorldToIndex(wc0, idx); mitk::Point3D cIdx; cIdx[0]=idx[0]; cIdx[1]=idx[1]; cIdx[2]=idx[2]; mitk::Point3D world; geom->IndexToWorld(cIdx,world); mitk::Vector3D trans = world - wc0; pe->GetGeometry()->Translate(trans); } } } } } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnFlipButton() { if (m_SelectedFiducial.IsNull()) return; std::map::iterator it = m_DataNodeToPlanarFigureData.find(m_SelectedFiducial.GetPointer()); if( it != m_DataNodeToPlanarFigureData.end() ) { QmitkPlanarFigureData& data = it->second; data.m_Flipped += 1; data.m_Flipped %= 2; } if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } QmitkFiberfoxView::GradientListType QmitkFiberfoxView::GenerateHalfShell(int NPoints) { NPoints *= 2; GradientListType pointshell; int numB0 = NPoints/20; if (numB0==0) numB0=1; GradientType g; g.Fill(0.0); for (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 > QmitkFiberfoxView::MakeGradientList() { std::vector > retval; vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); // Add 0 vector for B0 int numB0 = ndirs/10; if (numB0==0) numB0=1; itk::Vector v; v.Fill(0.0); for (int i=0; i v; v[0] = U->get(0,i); v[1] = U->get(1,i); v[2] = U->get(2,i); retval.push_back(v); } return retval; } void QmitkFiberfoxView::OnAddBundle() { if (m_SelectedImageNode.IsNull()) return; mitk::DataStorage::SetOfObjects::ConstPointer children = GetDataStorage()->GetDerivations(m_SelectedImageNode); mitk::FiberBundle::Pointer bundle = mitk::FiberBundle::New(); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( bundle ); QString name = QString("Bundle_%1").arg(children->size()); node->SetName(name.toStdString()); m_SelectedBundles.push_back(node); UpdateGui(); GetDataStorage()->Add(node, m_SelectedImageNode); } void QmitkFiberfoxView::OnDrawROI() { if (m_SelectedBundles.empty()) OnAddBundle(); if (m_SelectedBundles.empty()) return; mitk::DataStorage::SetOfObjects::ConstPointer children = GetDataStorage()->GetDerivations(m_SelectedBundles.at(0)); mitk::PlanarEllipse::Pointer figure = mitk::PlanarEllipse::New(); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( figure ); node->SetBoolProperty("planarfigure.3drendering", true); node->SetBoolProperty("planarfigure.3drendering.fill", true); QList nodes = this->GetDataManagerSelection(); for( int i=0; iSetSelected(false); m_SelectedFiducial = node; QString name = QString("Fiducial_%1").arg(children->size()); node->SetName(name.toStdString()); node->SetSelected(true); this->DisableCrosshairNavigation(); mitk::PlanarFigureInteractor::Pointer figureInteractor = dynamic_cast(node->GetDataInteractor().GetPointer()); if(figureInteractor.IsNull()) { figureInteractor = mitk::PlanarFigureInteractor::New(); us::Module* planarFigureModule = us::ModuleRegistry::GetModule( "MitkPlanarFigure" ); figureInteractor->LoadStateMachine("PlanarFigureInteraction.xml", planarFigureModule ); figureInteractor->SetEventConfig( "PlanarFigureConfig.xml", planarFigureModule ); figureInteractor->SetDataNode( node ); } UpdateGui(); GetDataStorage()->Add(node, m_SelectedBundles.at(0)); } bool CompareLayer(mitk::DataNode::Pointer i,mitk::DataNode::Pointer j) { int li = -1; i->GetPropertyValue("layer", li); int lj = -1; j->GetPropertyValue("layer", lj); return liGetSources(m_SelectedFiducial); for( mitk::DataStorage::SetOfObjects::const_iterator it = parents->begin(); it != parents->end(); ++it ) if(dynamic_cast((*it)->GetData())) m_SelectedBundles.push_back(*it); if (m_SelectedBundles.empty()) return; } FiberfoxParameters parameters; // = UpdateImageParameters(); for (unsigned int i=0; iGetDerivations(m_SelectedBundles.at(i)); std::vector< mitk::DataNode::Pointer > childVector; for( mitk::DataStorage::SetOfObjects::const_iterator it = children->begin(); it != children->end(); ++it ) childVector.push_back(*it); sort(childVector.begin(), childVector.end(), CompareLayer); vector< mitk::PlanarEllipse::Pointer > fib; vector< unsigned int > flip; float radius = 1; int count = 0; for( std::vector< mitk::DataNode::Pointer >::const_iterator it = childVector.begin(); it != childVector.end(); ++it ) { mitk::DataNode::Pointer node = *it; if ( node.IsNotNull() && dynamic_cast(node->GetData()) ) { mitk::PlanarEllipse* ellipse = dynamic_cast(node->GetData()); if (m_Controls->m_ConstantRadiusBox->isChecked()) { ellipse->SetTreatAsCircle(true); mitk::Point2D c = ellipse->GetControlPoint(0); mitk::Point2D p = ellipse->GetControlPoint(1); mitk::Vector2D v = p-c; if (count==0) { radius = v.GetVnlVector().magnitude(); ellipse->SetControlPoint(1, p); ellipse->Modified(); } else { v.Normalize(); v *= radius; ellipse->SetControlPoint(1, c+v); ellipse->Modified(); } } fib.push_back(ellipse); std::map::iterator it = m_DataNodeToPlanarFigureData.find(node.GetPointer()); if( it != m_DataNodeToPlanarFigureData.end() ) { QmitkPlanarFigureData& data = it->second; flip.push_back(data.m_Flipped); } else flip.push_back(0); } count++; } if (fib.size()>1) { parameters.m_FiberGen.m_Fiducials.push_back(fib); parameters.m_FiberGen.m_FlipList.push_back(flip); } else if (fib.size()>0) m_SelectedBundles.at(i)->SetData( mitk::FiberBundle::New() ); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } itk::FibersFromPlanarFiguresFilter::Pointer filter = itk::FibersFromPlanarFiguresFilter::New(); filter->SetParameters(parameters.m_FiberGen); filter->Update(); vector< mitk::FiberBundle::Pointer > fiberBundles = filter->GetFiberBundles(); for (unsigned int i=0; iSetData( fiberBundles.at(i) ); if (fiberBundles.at(i)->GetNumFibers()>50000) m_SelectedBundles.at(i)->SetVisibility(false); } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberfoxView::GenerateImage() { if (m_Controls->m_FiberBundleComboBox->GetSelectedNode().IsNull() && !mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( m_Controls->m_TemplateComboBox->GetSelectedNode())) { mitk::Image::Pointer image = mitk::ImageGenerator::GenerateGradientImage( m_Controls->m_SizeX->value(), m_Controls->m_SizeY->value(), m_Controls->m_SizeZ->value(), m_Controls->m_SpacingX->value(), m_Controls->m_SpacingY->value(), m_Controls->m_SpacingZ->value()); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); node->SetName("Dummy"); unsigned int window = m_Controls->m_SizeX->value()*m_Controls->m_SizeY->value()*m_Controls->m_SizeZ->value(); unsigned int level = window/2; mitk::LevelWindow lw; lw.SetLevelWindow(level, window); node->SetProperty( "levelwindow", mitk::LevelWindowProperty::New( lw ) ); GetDataStorage()->Add(node); m_SelectedImageNode = node; mitk::BaseData::Pointer basedata = node->GetData(); if (basedata.IsNotNull()) { mitk::RenderingManager::GetInstance()->InitializeViews( basedata->GetTimeGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true ); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } UpdateGui(); QMessageBox::information(NULL, "Template image generated", "You have selected no fiber bundle or diffusion-weighted image, which can be used to simulate a new diffusion-weighted image. A template image with the specified geometry has been generated that can be used to draw artificial fibers (see tab 'Fiber Definition')."); } else if (m_Controls->m_FiberBundleComboBox->GetSelectedNode().IsNotNull()) SimulateImageFromFibers(m_Controls->m_FiberBundleComboBox->GetSelectedNode()); else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( m_Controls->m_TemplateComboBox->GetSelectedNode()) ) SimulateForExistingDwi(m_Controls->m_TemplateComboBox->GetSelectedNode()); else QMessageBox::information(NULL, "No image generated", "You have selected no fiber bundle or diffusion-weighted image, which can be used to simulate a new diffusion-weighted image."); } void QmitkFiberfoxView::SimulateForExistingDwi(mitk::DataNode* imageNode) { bool isDiffusionImage( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(imageNode->GetData())) ); if ( !isDiffusionImage ) { return; } FiberfoxParameters parameters = UpdateImageParameters(); if (parameters.m_NoiseModel==NULL && parameters.m_SignalGen.m_Spikes==0 && parameters.m_SignalGen.m_FrequencyMap.IsNull() && parameters.m_SignalGen.m_KspaceLineOffset<=0.000001 && !parameters.m_SignalGen.m_DoAddGibbsRinging && !(parameters.m_SignalGen.m_EddyStrength>0) && parameters.m_SignalGen.m_CroppingFactor>0.999) { QMessageBox::information( NULL, "Simulation cancelled", "No valid artifact enabled! Motion artifacts and relaxation effects can NOT be added to an existing diffusion weighted image."); return; } mitk::Image::Pointer diffImg = dynamic_cast(imageNode->GetData()); ItkDwiType::Pointer itkVectorImagePointer = ItkDwiType::New(); mitk::CastToItkImage(diffImg, itkVectorImagePointer); m_ArtifactsToDwiFilter = itk::AddArtifactsToDwiImageFilter< short >::New(); m_ArtifactsToDwiFilter->SetInput(itkVectorImagePointer); parameters.m_Misc.m_ParentNode = imageNode; m_ArtifactsToDwiFilter->SetParameters(parameters); m_Worker.m_FilterType = 1; m_Thread.start(QThread::LowestPriority); } void QmitkFiberfoxView::SimulateImageFromFibers(mitk::DataNode* fiberNode) { mitk::FiberBundle::Pointer fiberBundle = dynamic_cast(fiberNode->GetData()); if (fiberBundle->GetNumFibers()<=0) return; FiberfoxParameters parameters = UpdateImageParameters(); m_TractsToDwiFilter = itk::TractsToDWIImageFilter< short >::New(); parameters.m_Misc.m_ParentNode = fiberNode; if (m_Controls->m_TemplateComboBox->GetSelectedNode().IsNotNull() && mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_Controls->m_TemplateComboBox->GetSelectedNode()->GetData()))) { bool first = true; bool ok = true; mitk::Image::Pointer diffImg = dynamic_cast(m_Controls->m_TemplateComboBox->GetSelectedNode()->GetData()); itk::Image< itk::DiffusionTensor3D< double >, 3 >::Pointer tensorImage = NULL; const int shOrder = 2; typedef itk::AnalyticalDiffusionQballReconstructionImageFilter QballFilterType; QballFilterType::CoefficientImageType::Pointer itkFeatureImage = NULL; ItkDoubleImgType::Pointer adcImage = NULL; for (unsigned int i=0; i* model = NULL; if (i* >(parameters.m_FiberModelList.at(i)); else model = dynamic_cast< mitk::RawShModel<>* >(parameters.m_NonFiberModelList.at(i-parameters.m_FiberModelList.size())); if (model!=0 && model->GetNumberOfKernels()<=0) { if (first==true) { ItkDwiType::Pointer itkVectorImagePointer = ItkDwiType::New(); mitk::CastToItkImage(diffImg, itkVectorImagePointer); typedef itk::DiffusionTensor3DReconstructionImageFilter< short, short, double > TensorReconstructionImageFilterType; TensorReconstructionImageFilterType::Pointer filter = TensorReconstructionImageFilterType::New(); filter->SetGradientImage( static_cast( diffImg->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer(), itkVectorImagePointer ); filter->SetBValue( static_cast(diffImg->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue() ); filter->Update(); tensorImage = filter->GetOutput(); const int NumCoeffs = (shOrder*shOrder + shOrder + 2)/2 + shOrder; QballFilterType::Pointer qballfilter = QballFilterType::New(); qballfilter->SetGradientImage( static_cast( diffImg->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer(), itkVectorImagePointer ); qballfilter->SetBValue( static_cast(diffImg->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue() ); qballfilter->SetLambda(0.006); qballfilter->SetNormalizationMethod(QballFilterType::QBAR_RAW_SIGNAL); qballfilter->Update(); itkFeatureImage = qballfilter->GetCoefficientImage(); itk::AdcImageFilter< short, double >::Pointer adcFilter = itk::AdcImageFilter< short, double >::New(); adcFilter->SetInput( itkVectorImagePointer ); adcFilter->SetGradientDirections( static_cast( diffImg->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer() ); adcFilter->SetB_value( static_cast(diffImg->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue() ); adcFilter->Update(); adcImage = adcFilter->GetOutput(); } ok = model->SampleKernels(diffImg, parameters.m_SignalGen.m_MaskImage, tensorImage, itkFeatureImage, adcImage); if (!ok) break; } } if (!ok) { QMessageBox::information( NULL, "Simulation cancelled", "No valid prototype signals could be sampled."); return; } } else if ( m_Controls->m_Compartment1Box->currentIndex()==3 || m_Controls->m_Compartment3Box->currentIndex()==3 || m_Controls->m_Compartment4Box->currentIndex()==4 ) { QMessageBox::information( NULL, "Simulation cancelled", "Prototype signal but no diffusion-weighted image selected to sample signal from."); return; } m_TractsToDwiFilter->SetParameters(parameters); m_TractsToDwiFilter->SetFiberBundle(fiberBundle); m_Worker.m_FilterType = 0; m_Thread.start(QThread::LowestPriority); } void QmitkFiberfoxView::ApplyTransform() { vector< mitk::DataNode::Pointer > selectedBundles; for(unsigned int i=0; iGetDerivations(m_SelectedImages.at(i)); for( mitk::DataStorage::SetOfObjects::const_iterator it = derivations->begin(); it != derivations->end(); ++it ) { mitk::DataNode::Pointer fibNode = *it; if ( fibNode.IsNotNull() && dynamic_cast(fibNode->GetData()) ) selectedBundles.push_back(fibNode); } } if (selectedBundles.empty()) selectedBundles = m_SelectedBundles2; if (!selectedBundles.empty()) { for (std::vector::const_iterator it = selectedBundles.begin(); it!=selectedBundles.end(); ++it) { mitk::FiberBundle::Pointer fib = dynamic_cast((*it)->GetData()); fib->RotateAroundAxis(m_Controls->m_XrotBox->value(), m_Controls->m_YrotBox->value(), m_Controls->m_ZrotBox->value()); fib->TranslateFibers(m_Controls->m_XtransBox->value(), m_Controls->m_YtransBox->value(), m_Controls->m_ZtransBox->value()); fib->ScaleFibers(m_Controls->m_XscaleBox->value(), m_Controls->m_YscaleBox->value(), m_Controls->m_ZscaleBox->value()); // handle child fiducials if (m_Controls->m_IncludeFiducials->isChecked()) { mitk::DataStorage::SetOfObjects::ConstPointer derivations = GetDataStorage()->GetDerivations(*it); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = derivations->begin(); it2 != derivations->end(); ++it2 ) { mitk::DataNode::Pointer fiducialNode = *it2; if ( fiducialNode.IsNotNull() && dynamic_cast(fiducialNode->GetData()) ) { mitk::PlanarEllipse* pe = dynamic_cast(fiducialNode->GetData()); mitk::BaseGeometry* geom = pe->GetGeometry(); // translate mitk::Vector3D world; world[0] = m_Controls->m_XtransBox->value(); world[1] = m_Controls->m_YtransBox->value(); world[2] = m_Controls->m_ZtransBox->value(); geom->Translate(world); // calculate rotation matrix double x = m_Controls->m_XrotBox->value()*M_PI/180; double y = m_Controls->m_YrotBox->value()*M_PI/180; double z = m_Controls->m_ZrotBox->value()*M_PI/180; itk::Matrix< double, 3, 3 > rotX; rotX.SetIdentity(); rotX[1][1] = cos(x); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(x); rotX[2][1] = -rotX[1][2]; itk::Matrix< double, 3, 3 > rotY; rotY.SetIdentity(); rotY[0][0] = cos(y); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(y); rotY[2][0] = -rotY[0][2]; itk::Matrix< double, 3, 3 > rotZ; rotZ.SetIdentity(); rotZ[0][0] = cos(z); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(z); rotZ[1][0] = -rotZ[0][1]; itk::Matrix< double, 3, 3 > rot = rotZ*rotY*rotX; // transform control point coordinate into geometry translation geom->SetOrigin(pe->GetWorldControlPoint(0)); mitk::Point2D cp; cp.Fill(0.0); pe->SetControlPoint(0, cp); // rotate fiducial geom->GetIndexToWorldTransform()->SetMatrix(rot*geom->GetIndexToWorldTransform()->GetMatrix()); // implicit translation mitk::Vector3D trans; trans[0] = geom->GetOrigin()[0]-fib->GetGeometry()->GetCenter()[0]; trans[1] = geom->GetOrigin()[1]-fib->GetGeometry()->GetCenter()[1]; trans[2] = geom->GetOrigin()[2]-fib->GetGeometry()->GetCenter()[2]; mitk::Vector3D newWc = rot*trans; newWc = newWc-trans; geom->Translate(newWc); pe->Modified(); } } } } } else { for (unsigned int i=0; i(m_SelectedFiducials.at(i)->GetData()); mitk::BaseGeometry* geom = pe->GetGeometry(); // translate mitk::Vector3D world; world[0] = m_Controls->m_XtransBox->value(); world[1] = m_Controls->m_YtransBox->value(); world[2] = m_Controls->m_ZtransBox->value(); geom->Translate(world); // calculate rotation matrix double x = m_Controls->m_XrotBox->value()*M_PI/180; double y = m_Controls->m_YrotBox->value()*M_PI/180; double z = m_Controls->m_ZrotBox->value()*M_PI/180; itk::Matrix< double, 3, 3 > rotX; rotX.SetIdentity(); rotX[1][1] = cos(x); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(x); rotX[2][1] = -rotX[1][2]; itk::Matrix< double, 3, 3 > rotY; rotY.SetIdentity(); rotY[0][0] = cos(y); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(y); rotY[2][0] = -rotY[0][2]; itk::Matrix< double, 3, 3 > rotZ; rotZ.SetIdentity(); rotZ[0][0] = cos(z); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(z); rotZ[1][0] = -rotZ[0][1]; itk::Matrix< double, 3, 3 > rot = rotZ*rotY*rotX; // transform control point coordinate into geometry translation geom->SetOrigin(pe->GetWorldControlPoint(0)); mitk::Point2D cp; cp.Fill(0.0); pe->SetControlPoint(0, cp); // rotate fiducial geom->GetIndexToWorldTransform()->SetMatrix(rot*geom->GetIndexToWorldTransform()->GetMatrix()); pe->Modified(); } if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberfoxView::CopyBundles() { if ( m_SelectedBundles.size()<1 ){ QMessageBox::information( NULL, "Warning", "Select at least one fiber bundle!"); MITK_WARN("QmitkFiberFoxView") << "Select at least one fiber bundle!"; return; } for (std::vector::const_iterator it = m_SelectedBundles.begin(); it!=m_SelectedBundles.end(); ++it) { // find parent image mitk::DataNode::Pointer parentNode; mitk::DataStorage::SetOfObjects::ConstPointer parentImgs = GetDataStorage()->GetSources(*it); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = parentImgs->begin(); it2 != parentImgs->end(); ++it2 ) { mitk::DataNode::Pointer pImgNode = *it2; if ( pImgNode.IsNotNull() && dynamic_cast(pImgNode->GetData()) ) { parentNode = pImgNode; break; } } mitk::FiberBundle::Pointer fib = dynamic_cast((*it)->GetData()); mitk::FiberBundle::Pointer newBundle = fib->GetDeepCopy(); QString name((*it)->GetName().c_str()); name += "_copy"; mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); if (parentNode.IsNotNull()) GetDataStorage()->Add(fbNode, parentNode); else GetDataStorage()->Add(fbNode); // copy child fiducials if (m_Controls->m_IncludeFiducials->isChecked()) { mitk::DataStorage::SetOfObjects::ConstPointer derivations = GetDataStorage()->GetDerivations(*it); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = derivations->begin(); it2 != derivations->end(); ++it2 ) { mitk::DataNode::Pointer fiducialNode = *it2; if ( fiducialNode.IsNotNull() && dynamic_cast(fiducialNode->GetData()) ) { mitk::PlanarEllipse::Pointer pe = dynamic_cast(fiducialNode->GetData())->Clone(); mitk::DataNode::Pointer newNode = mitk::DataNode::New(); newNode->SetData(pe); newNode->SetName(fiducialNode->GetName()); newNode->SetBoolProperty("planarfigure.3drendering", true); GetDataStorage()->Add(newNode, fbNode); } } } } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberfoxView::JoinBundles() { if ( m_SelectedBundles.size()<2 ){ QMessageBox::information( NULL, "Warning", "Select at least two fiber bundles!"); MITK_WARN("QmitkFiberFoxView") << "Select at least two fiber bundles!"; return; } std::vector::const_iterator it = m_SelectedBundles.begin(); mitk::FiberBundle::Pointer newBundle = dynamic_cast((*it)->GetData()); QString name(""); name += QString((*it)->GetName().c_str()); ++it; for (; it!=m_SelectedBundles.end(); ++it) { newBundle = newBundle->AddBundle(dynamic_cast((*it)->GetData())); name += "+"+QString((*it)->GetName().c_str()); } mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberfoxView::UpdateGui() { m_Controls->m_GeometryFrame->setEnabled(true); m_Controls->m_GeometryMessage->setVisible(false); m_Controls->m_DiffusionPropsMessage->setVisible(false); m_Controls->m_FiberGenMessage->setVisible(true); m_Controls->m_TransformBundlesButton->setEnabled(false); m_Controls->m_CopyBundlesButton->setEnabled(false); m_Controls->m_GenerateFibersButton->setEnabled(false); m_Controls->m_FlipButton->setEnabled(false); m_Controls->m_CircleButton->setEnabled(false); m_Controls->m_BvalueBox->setEnabled(true); m_Controls->m_NumGradientsBox->setEnabled(true); m_Controls->m_JoinBundlesButton->setEnabled(false); m_Controls->m_AlignOnGrid->setEnabled(false); // Fiber generation gui if (m_SelectedFiducial.IsNotNull()) { m_Controls->m_TransformBundlesButton->setEnabled(true); m_Controls->m_FlipButton->setEnabled(true); m_Controls->m_AlignOnGrid->setEnabled(true); } if (m_SelectedImageNode.IsNotNull() || !m_SelectedBundles.empty()) { m_Controls->m_CircleButton->setEnabled(true); m_Controls->m_FiberGenMessage->setVisible(false); } if (m_SelectedImageNode.IsNotNull() && !m_SelectedBundles.empty()) m_Controls->m_AlignOnGrid->setEnabled(true); if (!m_SelectedBundles.empty()) { m_Controls->m_TransformBundlesButton->setEnabled(true); m_Controls->m_CopyBundlesButton->setEnabled(true); m_Controls->m_GenerateFibersButton->setEnabled(true); if (m_SelectedBundles.size()>1) m_Controls->m_JoinBundlesButton->setEnabled(true); } // Signal generation gui if (m_Controls->m_MaskComboBox->GetSelectedNode().IsNotNull() || m_Controls->m_TemplateComboBox->GetSelectedNode().IsNotNull()) { m_Controls->m_GeometryMessage->setVisible(true); m_Controls->m_GeometryFrame->setEnabled(false); } if (m_Controls->m_TemplateComboBox->GetSelectedNode().IsNotNull() && mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_Controls->m_TemplateComboBox->GetSelectedNode()->GetData()))) { m_Controls->m_DiffusionPropsMessage->setVisible(true); m_Controls->m_BvalueBox->setEnabled(false); m_Controls->m_NumGradientsBox->setEnabled(false); m_Controls->m_GeometryMessage->setVisible(true); m_Controls->m_GeometryFrame->setEnabled(false); } } void QmitkFiberfoxView::OnSelectionChanged( berry::IWorkbenchPart::Pointer, const QList& nodes ) { m_SelectedBundles2.clear(); m_SelectedImages.clear(); m_SelectedFiducials.clear(); m_SelectedFiducial = NULL; m_SelectedBundles.clear(); m_SelectedImageNode = NULL; // iterate all selected objects, adjust warning visibility for( int i=0; i(node->GetData())); // } // if ( node.IsNotNull() && isDiffusionImage ) // { // m_SelectedDWI = node; // m_SelectedImage = node; // m_SelectedImages.push_back(node); // } if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { m_SelectedImages.push_back(node); m_SelectedImageNode = node; } else if ( node.IsNotNull() && dynamic_cast(node->GetData()) ) { m_SelectedBundles2.push_back(node); if (m_Controls->m_RealTimeFibers->isChecked()) { m_SelectedBundles.push_back(node); mitk::FiberBundle::Pointer newFib = dynamic_cast(node->GetData()); if (newFib->GetNumFibers()!=m_Controls->m_FiberDensityBox->value()) GenerateFibers(); } else m_SelectedBundles.push_back(node); } else if ( node.IsNotNull() && dynamic_cast(node->GetData()) ) { m_SelectedFiducials.push_back(node); m_SelectedFiducial = node; m_SelectedBundles.clear(); mitk::DataStorage::SetOfObjects::ConstPointer parents = GetDataStorage()->GetSources(node); for( mitk::DataStorage::SetOfObjects::const_iterator it = parents->begin(); it != parents->end(); ++it ) { mitk::DataNode::Pointer pNode = *it; if ( pNode.IsNotNull() && dynamic_cast(pNode->GetData()) ) m_SelectedBundles.push_back(pNode); } } } UpdateGui(); } void QmitkFiberfoxView::EnableCrosshairNavigation() { MITK_DEBUG << "EnableCrosshairNavigation"; // enable the crosshair navigation if (mitk::ILinkedRenderWindowPart* linkedRenderWindow = dynamic_cast(this->GetRenderWindowPart())) { MITK_DEBUG << "enabling linked navigation"; linkedRenderWindow->EnableLinkedNavigation(true); // linkedRenderWindow->EnableSlicingPlanes(true); } if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::DisableCrosshairNavigation() { MITK_DEBUG << "DisableCrosshairNavigation"; // disable the crosshair navigation during the drawing if (mitk::ILinkedRenderWindowPart* linkedRenderWindow = dynamic_cast(this->GetRenderWindowPart())) { MITK_DEBUG << "disabling linked navigation"; linkedRenderWindow->EnableLinkedNavigation(false); // linkedRenderWindow->EnableSlicingPlanes(false); } } void QmitkFiberfoxView::NodeRemoved(const mitk::DataNode* node) { mitk::DataNode* nonConstNode = const_cast(node); std::map::iterator it = m_DataNodeToPlanarFigureData.find(nonConstNode); if (dynamic_cast(node->GetData())) { m_SelectedBundles.clear(); m_SelectedBundles2.clear(); } else if (dynamic_cast(node->GetData())) m_SelectedImages.clear(); if( it != m_DataNodeToPlanarFigureData.end() ) { QmitkPlanarFigureData& data = it->second; // remove observers data.m_Figure->RemoveObserver( data.m_EndPlacementObserverTag ); data.m_Figure->RemoveObserver( data.m_SelectObserverTag ); data.m_Figure->RemoveObserver( data.m_StartInteractionObserverTag ); data.m_Figure->RemoveObserver( data.m_EndInteractionObserverTag ); m_DataNodeToPlanarFigureData.erase( it ); } } void QmitkFiberfoxView::NodeAdded( const mitk::DataNode* node ) { // add observer for selection in renderwindow mitk::PlanarFigure* figure = dynamic_cast(node->GetData()); bool isPositionMarker (false); node->GetBoolProperty("isContourMarker", isPositionMarker); if( figure && !isPositionMarker ) { MITK_DEBUG << "figure added. will add interactor if needed."; mitk::PlanarFigureInteractor::Pointer figureInteractor = dynamic_cast(node->GetDataInteractor().GetPointer()); mitk::DataNode* nonConstNode = const_cast( node ); if(figureInteractor.IsNull()) { figureInteractor = mitk::PlanarFigureInteractor::New(); us::Module* planarFigureModule = us::ModuleRegistry::GetModule( "MitkPlanarFigure" ); figureInteractor->LoadStateMachine("PlanarFigureInteraction.xml", planarFigureModule ); figureInteractor->SetEventConfig( "PlanarFigureConfig.xml", planarFigureModule ); figureInteractor->SetDataNode( nonConstNode ); } MITK_DEBUG << "will now add observers for planarfigure"; QmitkPlanarFigureData data; data.m_Figure = figure; // // add observer for event when figure has been placed typedef itk::SimpleMemberCommand< QmitkFiberfoxView > SimpleCommandType; // SimpleCommandType::Pointer initializationCommand = SimpleCommandType::New(); // initializationCommand->SetCallbackFunction( this, &QmitkFiberfoxView::PlanarFigureInitialized ); // data.m_EndPlacementObserverTag = figure->AddObserver( mitk::EndPlacementPlanarFigureEvent(), initializationCommand ); // add observer for event when figure is picked (selected) typedef itk::MemberCommand< QmitkFiberfoxView > MemberCommandType; MemberCommandType::Pointer selectCommand = MemberCommandType::New(); selectCommand->SetCallbackFunction( this, &QmitkFiberfoxView::PlanarFigureSelected ); data.m_SelectObserverTag = figure->AddObserver( mitk::SelectPlanarFigureEvent(), selectCommand ); // add observer for event when interaction with figure starts SimpleCommandType::Pointer startInteractionCommand = SimpleCommandType::New(); startInteractionCommand->SetCallbackFunction( this, &QmitkFiberfoxView::DisableCrosshairNavigation); data.m_StartInteractionObserverTag = figure->AddObserver( mitk::StartInteractionPlanarFigureEvent(), startInteractionCommand ); // add observer for event when interaction with figure starts SimpleCommandType::Pointer endInteractionCommand = SimpleCommandType::New(); endInteractionCommand->SetCallbackFunction( this, &QmitkFiberfoxView::EnableCrosshairNavigation); data.m_EndInteractionObserverTag = figure->AddObserver( mitk::EndInteractionPlanarFigureEvent(), endInteractionCommand ); m_DataNodeToPlanarFigureData[nonConstNode] = data; } } void QmitkFiberfoxView::PlanarFigureSelected( itk::Object* object, const itk::EventObject& ) { mitk::TNodePredicateDataType::Pointer isPf = mitk::TNodePredicateDataType::New(); mitk::DataStorage::SetOfObjects::ConstPointer allPfs = this->GetDataStorage()->GetSubset( isPf ); for ( mitk::DataStorage::SetOfObjects::const_iterator it = allPfs->begin(); it!=allPfs->end(); ++it) { mitk::DataNode* node = *it; if( node->GetData() == object ) { node->SetSelected(true); m_SelectedFiducial = node; } else node->SetSelected(false); } UpdateGui(); this->RequestRenderWindowUpdate(); } void QmitkFiberfoxView::SetFocus() { m_Controls->m_CircleButton->setFocus(); } void QmitkFiberfoxView::SetOutputPath() { // SELECT FOLDER DIALOG string outputPath = QFileDialog::getExistingDirectory(NULL, "Save images to...", QString(outputPath.c_str())).toStdString(); if (outputPath.empty()) m_Controls->m_SavePathEdit->setText("-"); else { outputPath += "/"; m_Controls->m_SavePathEdit->setText(QString(outputPath.c_str())); } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxViewControls.ui index 6afa25f654..b1cb739ca5 100755 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxViewControls.ui @@ -1,3254 +1,3295 @@ QmitkFiberfoxViewControls 0 0 524 2652 Form Load Parameters :/QmitkDiffusionImaging/general_icons/upload.ico:/QmitkDiffusionImaging/general_icons/upload.ico - 1 + 0 Fiber Definition Qt::Vertical 20 40 color: rgb(255, 0, 0); Please select an image or an existing fiber bundle to draw the fiber fiducials. If you can't provide a suitable image, generate one using the "Signal Generation" tab. Qt::AutoText Qt::AlignJustify|Qt::AlignVCenter true Fiducial Options All fiducials are treated as circles with the same radius as the first fiducial. Use Constant Fiducial Radius false false Align selected fiducials with voxel grid. Shifts selected fiducials to nearest voxel center. Align With Grid :/QmitkDiffusionImaging/general_icons/right.ico:/QmitkDiffusionImaging/general_icons/right.ico Operations false Join Bundles :/QmitkDiffusionImaging/general_icons/plus.ico:/QmitkDiffusionImaging/general_icons/plus.ico QFrame::NoFrame QFrame::Raised 0 0 0 0 Y false Rotation angle (in degree) around x-axis. 3 -360.000000000000000 360.000000000000000 0.100000000000000 Axis: false Rotation angle (in degree) around y-axis. 3 -360.000000000000000 360.000000000000000 0.100000000000000 Translation: false Translation (in mm) in direction of the z-axis. 3 -1000.000000000000000 1000.000000000000000 0.100000000000000 Translation (in mm) in direction of the y-axis. 3 -1000.000000000000000 1000.000000000000000 0.100000000000000 X false Rotation: false Z false Rotation angle (in degree) around z-axis. 3 -360.000000000000000 360.000000000000000 0.100000000000000 Translation (in mm) in direction of the x-axis. 3 -1000.000000000000000 1000.000000000000000 0.100000000000000 Scaling: false Scaling factor for selected fiber bundle along the x-axis. 0.010000000000000 10.000000000000000 0.010000000000000 1.000000000000000 Scaling factor for selected fiber bundle along the y-axis. 0.010000000000000 10.000000000000000 0.010000000000000 1.000000000000000 Scaling factor for selected fiber bundle along the z-axis. 0.010000000000000 10.000000000000000 0.010000000000000 1.000000000000000 false Copy Bundles :/QmitkDiffusionImaging/general_icons/copy2.ico:/QmitkDiffusionImaging/general_icons/copy2.ico false Transform Selection :/QmitkDiffusionImaging/general_icons/refresh.ico:/QmitkDiffusionImaging/general_icons/refresh.ico If checked, the fiducials belonging to the modified bundle are also modified. Include Fiducials true Fiber Options QFrame::NoFrame QFrame::Raised 0 0 0 0 QFrame::NoFrame QFrame::Raised 0 0 0 0 Tension: false Fiber Sampling: false 3 -1.000000000000000 1.000000000000000 0.100000000000000 0.000000000000000 3 -1.000000000000000 1.000000000000000 0.100000000000000 0.000000000000000 Bias: false Continuity: false 3 -1.000000000000000 1.000000000000000 0.100000000000000 0.000000000000000 Distance of fiber sampling points (in mm) 1 0.100000000000000 0.100000000000000 1.000000000000000 QFrame::NoFrame QFrame::Raised 0 0 0 0 6 #Fibers: false Specify number of fibers to generate for the selected bundle. 1 1000000 100 100 false Generate Fibers :/QmitkDiffusionImaging/general_icons/right.ico:/QmitkDiffusionImaging/general_icons/right.ico QFrame::NoFrame QFrame::Raised 0 0 0 0 Select fiber distribution inside of the fiducials. Uniform Gaussian Fiber Distribution: false Variance of the gaussian 3 0.001000000000000 10.000000000000000 0.010000000000000 0.100000000000000 QFrame::NoFrame QFrame::Raised 0 0 0 0 Disable to only generate fibers if "Generate Fibers" button is pressed. Real Time Fibers true Disable to only generate fibers if "Generate Fibers" button is pressed. Advanced Options false QFrame::NoFrame QFrame::Raised 0 0 0 0 false 30 30 Draw elliptical fiducial. :/QmitkDiffusionImaging/circle.png:/QmitkDiffusionImaging/circle.png 32 32 false true false 30 30 Flip fiber waypoints of selcted fiducial around one axis. :/QmitkDiffusionImaging/refresh.xpm:/QmitkDiffusionImaging/refresh.xpm 32 32 false true Qt::Horizontal 40 20 Signal Generation Extra-axonal Compartments QFrame::NoFrame QFrame::Raised 0 0 0 0 Volume Fraction: Select signal model for extra-axonal compartment. Ball Model Astrosticks Model Dot Model Prototype Signal Qt::Horizontal QFrame::NoFrame QFrame::Raised 0 0 0 0 Volume Fraction: Optional! If no volume fraction map for this compartment is set, the corresponding volume fractions are calculated from the input fibers. Select signal model for extra-axonal compartment. -- Ball Model Astrosticks Model Dot Model Prototype Signal Qt::Vertical 20 40 Image Settings Advanced Options QFrame::NoFrame QFrame::Raised 0 0 0 0 6 Gradient Directions: Number of gradient directions distributed over the half sphere. 0 10000 1 30 <html><head/><body><p>b-Value<span style=" font-style:italic;"> [s/mm</span><span style=" font-style:italic; vertical-align:super;">2</span><span style=" font-style:italic;">]</span>:</p></body></html> false b-value in s/mm² 0 10000 100 1000 color: rgb(255, 0, 0); Using geometry of selected image! color: rgb(255, 0, 0); Using gradients of selected DWI! QFrame::NoFrame QFrame::Raised 0 0 0 0 6 - - + + - Output one image per compartment containing the corresponding volume fractions per voxel. + Disable partial volume. Treat voxel content as fiber-only if at least one fiber is present. - Reverse Phase Encoding Direction + Disable Partial Volume Effects false - - + + + + + + + + + + + - Signal Scale: + <html><head/><body><p>Repetition Time <span style=" font-style:italic;">TR</span>: </p></body></html> + + + false - - + + - Relaxation time due to magnetic field inhomogeneities (T2', in milliseconds). + TE in milliseconds 1 10000 1 - 50 - - - - - - - Output one image per compartment containing the corresponding volume fractions per voxel. - - - Output Volume Fractions - - - false + 100 Fiber Radius: + + + + Relaxation time due to magnetic field inhomogeneities (T2', in milliseconds). + + + 1 + + + 10000 + + + 1 + + + 50 + + + T2* relaxation time (in milliseconds). 100.000000000000000 0.100000000000000 1.000000000000000 - - + + - <html><head/><body><p><span style=" font-style:italic;">TE</span>, <span style=" font-style:italic;">T</span><span style=" font-style:italic; vertical-align:sub;">inhom</span> and <span style=" font-style:italic;">T2</span> will have no effect if unchecked.</p></body></html> + + + + + + + - Simulate Signal Relaxation + Partial Fourier: - - true + + false - - + + - Fiber radius used to calculate volume fractions (in µm). Set to 0 for automatic radius estimation. + + + + + + + + + + <html><head/><body><p>Echo Time <span style=" font-style:italic;">TE</span>: </p></body></html> + + + false + + + + + + + TR in milliseconds - 0 + 1 - 1000 + 999999999 + + + 1 - 0 + 2000 - - + + TE in milliseconds 1 10000 1 100 + + + + Output phase image and volume fraction maps. + + + Output Additional Images + + + false + + + Dwell time: false - - + + - Disable partial volume. Treat voxel content as fiber-only if at least one fiber is present. + Fiber radius used to calculate volume fractions (in µm). Set to 0 for automatic radius estimation. - - Disable Partial Volume Effects + + 0 - - false + + 1000 + + + 0 - - + + - - - - - - - + <html><head/><body><p><span style=" font-style:italic;">TE</span>, <span style=" font-style:italic;">T</span><span style=" font-style:italic; vertical-align:sub;">inhom</span> and <span style=" font-style:italic;">T2</span> will have no effect if unchecked.</p></body></html> - <html><head/><body><p>Echo Time <span style=" font-style:italic;">TE</span>: </p></body></html> + Simulate Signal Relaxation - - false + + true - - + + - TE in milliseconds + T2* relaxation time (in milliseconds). + + + 3 - 1 + 0.500000000000000 - 10000 + 1.000000000000000 - 1 + 0.100000000000000 - 100 + 1.000000000000000 <html><head/><body><p><span style=" font-style:italic;">T</span><span style=" font-style:italic; vertical-align:sub;">inhom</span> Relaxation: </p></body></html> false - - - - TR in milliseconds - - - 1 - - - 999999999 - - - 1 - - - 2000 + + + + Signal Scale: - - + + - - - - - - - + Output one image per compartment containing the corresponding volume fractions per voxel. - <html><head/><body><p>Repetition Time <span style=" font-style:italic;">TR</span>: </p></body></html> + Reverse Phase Encoding Direction - + false QFrame::NoFrame QFrame::Raised 0 0 0 0 3 0.100000000000000 50.000000000000000 0.100000000000000 2.000000000000000 Image Spacing: 3 0.100000000000000 50.000000000000000 0.100000000000000 2.000000000000000 3 0.100000000000000 50.000000000000000 0.100000000000000 2.000000000000000 Image Dimensions: Fiber sampling factor which determines the accuracy of the calculated fiber and non-fiber volume fractions. 1 1000 1 11 Fiber sampling factor which determines the accuracy of the calculated fiber and non-fiber volume fractions. 1 1000 1 11 Fiber sampling factor which determines the accuracy of the calculated fiber and non-fiber volume fractions. 1 1000 1 3 Inter-axonal Compartment Select signal model for intra-axonal compartment. -- Stick Model Zeppelin Model Tensor Model QFrame::NoFrame QFrame::Raised 0 0 0 0 Volume Fraction: Optional! If no volume fraction map for this compartment is set, the corresponding volume fractions are calculated from the input fibers. true Stop current simulation. Abort Simulation :/QmitkDiffusionImaging/general_icons/abort.ico:/QmitkDiffusionImaging/general_icons/abort.ico Data QFrame::NoFrame QFrame::Raised 0 0 0 0 0 - ... <html><head/><body><p>Select a binary image to define the area of signal generation. Outside of the mask image only noise will be actively generated.</p></body></html> QComboBox::AdjustToMinimumContentsLength Fiber Bundle: false Save path: false Tissue Mask: false <html><head/><body><p>Select a fiber bundle to generate the white matter signal from. You can either use the fiber definition tab to manually define an input fiber bundle or you can also use any existing bundle, e.g. yielded by a tractography algorithm.</p></body></html> QComboBox::AdjustToMinimumContentsLength Template Image: false <html><head/><body><p>The parameters for the simulation (e.g. spacing, size, diffuison-weighted gradients, b-value) are adopted from this image.</p></body></html> QComboBox::AdjustToMinimumContentsLength 8 true Noise and other Artifacts Qt::Horizontal Add Noise false Add ringing artifacts occuring at strong edges in the image. Add Gibbs Ringing false true QFrame::NoFrame QFrame::Raised 6 0 0 0 0 Shrink FOV (%): false Shrink FOV by this percentage. 1 0.000000000000000 90.000000000000000 0.100000000000000 25.000000000000000 Qt::Horizontal QFrame::NoFrame QFrame::Raised 0 0 0 0 Num. Spikes: The number of randomly occurring signal spikes. 1 Spike amplitude relative to the largest signal amplitude of the corresponding k-space slice. 0.100000000000000 0.100000000000000 Scale: !!!EXPERIMENTAL!!! Add Eddy Current Effects false Add Spikes false QFrame::NoFrame QFrame::Raised 0 0 0 0 Variance: Variance of selected noise distribution. - 4 + 10 0.000000000000000 999999999.000000000000000 0.001000000000000 50.000000000000000 Distribution: Noise distribution Rician Chi-squared Add N/2 Ghosts false true QFrame::NoFrame QFrame::Raised 6 0 0 0 0 Frequency Map: false Select image specifying the frequency inhomogeneities (in Hz). Qt::Horizontal Qt::Horizontal Qt::Horizontal true QFrame::NoFrame QFrame::Raised QFormLayout::AllNonFixedFieldsGrow 6 0 6 0 0 Toggle between random movement and linear movement. Randomize motion true Rotation 0 9 0 0 Degree: false x false Axis: false Maximum rotation around x-axis. 1 360.000000000000000 1.000000000000000 0.000000000000000 Maximum rotation around z-axis. 1 360.000000000000000 1.000000000000000 15.000000000000000 y false z false Maximum rotation around y-axis. 1 360.000000000000000 1.000000000000000 0.000000000000000 Translation 0 0 0 Distance: false x false y false Axis: false z false Maximum translation along x-axis. 1 1000.000000000000000 1.000000000000000 0.000000000000000 Maximum translation along y-axis. 1 1000.000000000000000 1.000000000000000 0.000000000000000 Maximum translation along z-axis. 1 1000.000000000000000 1.000000000000000 0.000000000000000 Add Motion Artifacts false Add Distortions false Add Aliasing false true QFrame::NoFrame QFrame::Raised 6 0 0 0 0 K-Space Line Offset: false A larger offset increases the inensity of the ghost image. 3 1.000000000000000 0.010000000000000 0.250000000000000 true QFrame::NoFrame QFrame::Raised QFormLayout::AllNonFixedFieldsGrow 6 0 0 0 0 Magnitude: false Maximum magnitude of eddy current induced magnetic field inhomogeneities (in mT). 5 1000.000000000000000 0.001000000000000 0.005000000000000 color: rgb(255, 0, 0); Experimental! Qt::Horizontal Qt::Horizontal true <html><head/><body><p>Start DWI generation from selected fiber bundle.</p><p>If no fiber bundle but an existing diffusion weighted image is selected, the enabled artifacts are added to this image.</p><p>If neither a fiber bundle nor a diffusion weighted image is selected, a grayscale image containing a simple gradient is generated.</p></body></html> Start Simulation :/QmitkDiffusionImaging/general_icons/right.ico:/QmitkDiffusionImaging/general_icons/right.ico Intra-axonal Compartment Select signal model for intra-axonal compartment. Stick Model Zeppelin Model Tensor Model Prototype Signal QFrame::NoFrame QFrame::Raised 0 0 0 0 Volume Fraction: Optional! If no volume fraction map for this compartment is set, the corresponding volume fractions are calculated from the input fibers. Save Parameters :/QmitkDiffusionImaging/general_icons/download.ico:/QmitkDiffusionImaging/general_icons/download.ico QmitkDataStorageComboBox QComboBox
QmitkDataStorageComboBox.h
 QmitkTensorModelParametersWidget QWidget
QmitkTensorModelParametersWidget.h
1 QmitkStickModelParametersWidget QWidget
QmitkStickModelParametersWidget.h
1 QmitkZeppelinModelParametersWidget QWidget
QmitkZeppelinModelParametersWidget.h
1 QmitkBallModelParametersWidget QWidget
QmitkBallModelParametersWidget.h
1 QmitkAstrosticksModelParametersWidget QWidget
QmitkAstrosticksModelParametersWidget.h
1 QmitkDotModelParametersWidget QWidget
QmitkDotModelParametersWidget.h
1 QmitkPrototypeSignalParametersWidget QWidget
QmitkPrototypeSignalParametersWidget.h
1 QmitkDataStorageComboBoxWithSelectNone QComboBox
QmitkDataStorageComboBoxWithSelectNone.h
 m_CircleButton m_FlipButton m_RealTimeFibers m_AdvancedOptionsBox m_DistributionBox m_VarianceBox m_FiberDensityBox m_FiberSamplingBox m_TensionBox m_ContinuityBox m_BiasBox m_GenerateFibersButton m_ConstantRadiusBox m_AlignOnGrid m_XrotBox m_YrotBox m_ZrotBox m_XtransBox m_YtransBox m_ZtransBox m_XscaleBox m_YscaleBox m_ZscaleBox m_TransformBundlesButton m_CopyBundlesButton m_JoinBundlesButton m_IncludeFiducials m_GenerateImageButton m_SizeX m_SizeY m_SizeZ m_SpacingX m_SpacingY m_SpacingZ m_NumGradientsBox m_BvalueBox m_AdvancedOptionsBox_2 m_SignalScaleBox m_TEbox m_LineReadoutTimeBox m_T2starBox m_FiberRadius m_RelaxationBox m_EnforcePureFiberVoxelsBox m_VolumeFractionsBox m_Compartment1Box m_Compartment2Box m_Compartment3Box m_Compartment4Box m_AddNoise m_NoiseLevel m_AddSpikes m_SpikeNumBox m_SpikeScaleBox m_AddGhosts m_kOffsetBox m_AddAliasing m_WrapBox m_AddDistortions m_FrequencyMapBox m_AddMotion m_RandomMotion m_MaxRotationBoxX m_MaxRotationBoxY m_MaxRotationBoxZ m_MaxTranslationBoxX m_MaxTranslationBoxY m_MaxTranslationBoxZ m_AddEddy m_EddyGradientStrength m_AddGibbsRinging m_SaveParametersButton m_LoadParametersButton tabWidget