diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/Sequences/mitkCartesianReadout.h b/Modules/DiffusionImaging/FiberTracking/Fiberfox/Sequences/mitkConventionalSpinEcho.h similarity index 85% rename from Modules/DiffusionImaging/FiberTracking/Fiberfox/Sequences/mitkCartesianReadout.h rename to Modules/DiffusionImaging/FiberTracking/Fiberfox/Sequences/mitkConventionalSpinEcho.h index 17073912d9..aa262f8c2d 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/Sequences/mitkCartesianReadout.h +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/Sequences/mitkConventionalSpinEcho.h @@ -1,89 +1,89 @@ /*=================================================================== 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_CartesianReadout_H -#define _MITK_CartesianReadout_H +#ifndef _MITK_ConventionalSpinEcho_H +#define _MITK_ConventionalSpinEcho_H #include namespace mitk { /** - * \brief + * \brief Conventional spin echo sequence. Cartesian readout. One echo and one excitation per k-space line. * */ -class CartesianReadout : public AcquisitionType +class ConventionalSpinEcho : public AcquisitionType { public: - CartesianReadout(FiberfoxParameters* parameters) : AcquisitionType(parameters) + ConventionalSpinEcho(FiberfoxParameters* parameters) : AcquisitionType(parameters) { kxMax = m_Parameters->m_SignalGen.m_CroppedRegion.GetSize(0); kyMax = m_Parameters->m_SignalGen.m_CroppedRegion.GetSize(1); dt = m_Parameters->m_SignalGen.m_tLine/kxMax; // time to read one k-space voxel // maximum echo at center of each line m_NegTEhalf = -dt*(kxMax-(int)kxMax%2)/2; } - ~CartesianReadout() override + ~ConventionalSpinEcho() override {} float GetTimeFromMaxEcho(itk::Index< 2 > index) override { float t = 0; t = m_NegTEhalf + (float)index[0]*dt; return t; } float GetRedoutTime(itk::Index< 2 > index) override { float t = 0; t = (float)index[0]*dt; return t; } itk::Index< 2 > GetActualKspaceIndex(itk::Index< 2 > index) override { // reverse phase if (!m_Parameters->m_SignalGen.m_ReversePhase) index[1] = kyMax-1-index[1]; return index; } void AdjustEchoTime() override { if ( m_Parameters->m_SignalGen.m_tEcho/2 < m_Parameters->m_SignalGen.m_tLine/2 ) { m_Parameters->m_SignalGen.m_tEcho = m_Parameters->m_SignalGen.m_tLine; MITK_WARN << "Echo time is too short! Time not sufficient to read slice. Automatically adjusted to " << m_Parameters->m_SignalGen.m_tEcho << " ms"; m_Parameters->m_Misc.m_AfterSimulationMessage += "Echo time was chosen too short! Time not sufficient to read slice. Internally adjusted to " + boost::lexical_cast(m_Parameters->m_SignalGen.m_tEcho) + " ms\n"; } } protected: float dt; int kxMax; int kyMax; }; } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkKspaceImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkKspaceImageFilter.cpp index 335241aba6..acaba8ddd3 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkKspaceImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkKspaceImageFilter.cpp @@ -1,463 +1,463 @@ /*=================================================================== 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 //#endif #include #include #include #include #include "itkKspaceImageFilter.h" #include #include #include #include -#include +#include #include namespace itk { template< class ScalarType > KspaceImageFilter< ScalarType >::KspaceImageFilter() : m_Z(0) , m_RandSeed(-1) , m_SpikesPerSlice(0) , m_IsBaseline(true) { m_DiffusionGradientDirection.Fill(0.0); m_CoilPosition.Fill(0.0); } template< class ScalarType > void KspaceImageFilter< ScalarType > ::BeforeThreadedGenerateData() { m_Spike = vcl_complex(0,0); m_SpikeLog = ""; m_TransX = -m_Translation[0]; m_TransY = -m_Translation[1]; m_TransZ = -m_Translation[2]; kxMax = m_Parameters->m_SignalGen.m_CroppedRegion.GetSize(0); kyMax = m_Parameters->m_SignalGen.m_CroppedRegion.GetSize(1); xMax = m_CompartmentImages.at(0)->GetLargestPossibleRegion().GetSize(0); // scanner coverage in x-direction yMax = m_CompartmentImages.at(0)->GetLargestPossibleRegion().GetSize(1); // scanner coverage in y-direction yMaxFov = yMax; if (m_Parameters->m_Misc.m_DoAddAliasing) yMaxFov *= m_Parameters->m_SignalGen.m_CroppingFactor; // actual FOV in y-direction (in x-direction FOV=xMax) yMaxFov_half = yMaxFov/2; numPix = kxMax*kyMax; float ringing_factor = static_cast(m_Parameters->m_SignalGen.m_ZeroRinging)/100.0; ringing_lines_x = static_cast(ceil(kxMax/2 * ringing_factor)); ringing_lines_y = static_cast(ceil(kyMax/2 * ringing_factor)); // Adjust noise variance since it is the intended variance in physical space and not in k-space: float noiseVar = m_Parameters->m_SignalGen.m_PartialFourier*m_Parameters->m_SignalGen.m_NoiseVariance/(kyMax*kxMax); m_RandGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); if (m_RandSeed>=0) // always generate the same random numbers? m_RandGen->SetSeed(m_RandSeed); else m_RandGen->SetSeed(); typename OutputImageType::Pointer outputImage = OutputImageType::New(); itk::ImageRegion<2> region; region.SetSize(0, m_Parameters->m_SignalGen.m_CroppedRegion.GetSize(0)); region.SetSize(1, m_Parameters->m_SignalGen.m_CroppedRegion.GetSize(1)); outputImage->SetLargestPossibleRegion( region ); outputImage->SetBufferedRegion( region ); outputImage->SetRequestedRegion( region ); outputImage->Allocate(); vcl_complex zero = vcl_complex(0, 0); outputImage->FillBuffer(zero); if (m_Parameters->m_SignalGen.m_NoiseVariance>0 && m_Parameters->m_Misc.m_DoAddNoise) { ImageRegionIterator< OutputImageType > oit(outputImage, outputImage->GetLargestPossibleRegion()); while( !oit.IsAtEnd() ) { oit.Set(vcl_complex(m_RandGen->GetNormalVariate(0, noiseVar), m_RandGen->GetNormalVariate(0, noiseVar))); ++oit; } } m_KSpaceImage = InputImageType::New(); m_KSpaceImage->SetLargestPossibleRegion( region ); m_KSpaceImage->SetBufferedRegion( region ); m_KSpaceImage->SetRequestedRegion( region ); m_KSpaceImage->Allocate(); m_KSpaceImage->FillBuffer(0.0); m_Gamma = 42576000; // Gyromagnetic ratio in Hz/T (1.5T) 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 * m_Gamma; m_IsBaseline = false; } this->SetNthOutput(0, outputImage); for (int i=0; i<3; i++) for (int j=0; j<3; j++) m_Transform[i][j] = m_Parameters->m_SignalGen.m_ImageDirection[i][j] * m_Parameters->m_SignalGen.m_ImageSpacing[j]/1000; float a = m_Parameters->m_SignalGen.m_ImageRegion.GetSize(0)*m_Parameters->m_SignalGen.m_ImageSpacing[0]; float b = m_Parameters->m_SignalGen.m_ImageRegion.GetSize(1)*m_Parameters->m_SignalGen.m_ImageSpacing[1]; float diagonal = sqrt(a*a+b*b)/1000; // image diagonal in m switch (m_Parameters->m_SignalGen.m_CoilSensitivityProfile) { case SignalGenerationParameters::COIL_CONSTANT: { m_CoilSensitivityFactor = 1; // same signal everywhere break; } case SignalGenerationParameters::COIL_LINEAR: { m_CoilSensitivityFactor = -1/diagonal; // about 50% of the signal in the image center remaining break; } case SignalGenerationParameters::COIL_EXPONENTIAL: { m_CoilSensitivityFactor = -log(0.1)/diagonal; // about 32% of the signal in the image center remaining break; } } switch (m_Parameters->m_SignalGen.m_AcquisitionType) { case SignalGenerationParameters::SingleShotEpi: m_ReadoutScheme = new mitk::SingleShotEpi(m_Parameters); break; - case SignalGenerationParameters::SpinEcho: - m_ReadoutScheme = new mitk::CartesianReadout(m_Parameters); + case SignalGenerationParameters::ConventionalSpinEcho: + m_ReadoutScheme = new mitk::ConventionalSpinEcho(m_Parameters); break; default: m_ReadoutScheme = new mitk::SingleShotEpi(m_Parameters); } m_ReadoutScheme->AdjustEchoTime(); m_MovedFmap = nullptr; if (m_Parameters->m_Misc.m_DoAddDistortions && m_Parameters->m_SignalGen.m_FrequencyMap.IsNotNull() && m_Parameters->m_SignalGen.m_DoAddMotion) { // we have to account for the head motion since this also moves our frequency map itk::LinearInterpolateImageFunction< itk::Image< float, 3 >, float >::Pointer fmapInterpolator; fmapInterpolator = itk::LinearInterpolateImageFunction< itk::Image< float, 3 >, float >::New(); fmapInterpolator->SetInputImage(m_Parameters->m_SignalGen.m_FrequencyMap); m_MovedFmap = itk::Image< ScalarType, 2 >::New(); m_MovedFmap->SetLargestPossibleRegion( m_CompartmentImages.at(0)->GetLargestPossibleRegion() ); m_MovedFmap->SetBufferedRegion( m_CompartmentImages.at(0)->GetLargestPossibleRegion() ); m_MovedFmap->SetRequestedRegion( m_CompartmentImages.at(0)->GetLargestPossibleRegion() ); m_MovedFmap->Allocate(); m_MovedFmap->FillBuffer(0); ImageRegionIterator< InputImageType > it(m_MovedFmap, m_MovedFmap->GetLargestPossibleRegion() ); while( !it.IsAtEnd() ) { itk::Image::IndexType index; index[0] = it.GetIndex()[0]; index[1] = it.GetIndex()[1]; index[2] = m_Zidx; itk::Point point3D; m_Parameters->m_SignalGen.m_FrequencyMap->TransformIndexToPhysicalPoint(index, point3D); m_FiberBundle->TransformPoint( point3D, m_RotationMatrix, m_TransX, m_TransY, m_TransZ ); it.Set(mitk::imv::GetImageValue(point3D, true, fmapInterpolator)); ++it; } } // calculate T1 relaxation (independent of actual readout) m_T1Relax.clear(); if ( m_Parameters->m_SignalGen.m_DoSimulateRelaxation) for (unsigned int i=0; im_SignalGen.m_tRep/m_T1[i])); // account for inversion pulse and TI if (m_Parameters->m_SignalGen.m_tInv > 0) relaxation *= (1.0-std::exp(std::log(2) - m_Parameters->m_SignalGen.m_tInv/m_T1[i])); m_T1Relax.push_back(relaxation); } } template< class ScalarType > float KspaceImageFilter< ScalarType >::CoilSensitivity(VectorType& pos) { // ************************************************************************* // Coil ring is moving with excited slice (FIX THIS SOMETIME) m_CoilPosition[2] = pos[2]; // ************************************************************************* switch (m_Parameters->m_SignalGen.m_CoilSensitivityProfile) { case SignalGenerationParameters::COIL_CONSTANT: return 1; case SignalGenerationParameters::COIL_LINEAR: { VectorType diff = pos-m_CoilPosition; float sens = diff.GetNorm()*m_CoilSensitivityFactor + 1; if (sens<0) sens = 0; return sens; } case SignalGenerationParameters::COIL_EXPONENTIAL: { VectorType diff = pos-m_CoilPosition; float dist = static_cast(diff.GetNorm()); return std::exp(-dist*m_CoilSensitivityFactor); } default: return 1; } } template< class ScalarType > void KspaceImageFilter< ScalarType > ::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; vcl_complex zero = vcl_complex(0, 0); while( !oit.IsAtEnd() ) { typename OutputImageType::IndexType out_idx = oit.GetIndex(); // get current k-space index (depends on the chosen k-space readout scheme) itk::Index< 2 > kIdx = m_ReadoutScheme->GetActualKspaceIndex(out_idx); // partial fourier if (kIdx[1]>kyMax*m_Parameters->m_SignalGen.m_PartialFourier) { outputImage->SetPixel(kIdx, zero); ++oit; continue; } // gibbs ringing by setting high frequencies to zero (alternative to using smaller k-space than input image space) if (m_Parameters->m_SignalGen.m_DoAddGibbsRinging && m_Parameters->m_SignalGen.m_ZeroRinging>0) { if (kIdx[0] < ringing_lines_x || kIdx[1] < ringing_lines_y || kIdx[0] >= kxMax - ringing_lines_x || kIdx[1] >= kyMax - ringing_lines_y) { outputImage->SetPixel(kIdx, zero); ++oit; continue; } } // shift k for DFT: (0 -- N) --> (-N/2 -- N/2) float kx = kIdx[0]; float ky = kIdx[1]; if (static_cast(kxMax)%2==1) kx -= (kxMax-1)/2; else kx -= kxMax/2; if (static_cast(kyMax)%2==1) ky -= (kyMax-1)/2; else ky -= kyMax/2; // time from maximum echo float t = m_ReadoutScheme->GetTimeFromMaxEcho(out_idx); // time passed since k-space readout started float tRead = m_ReadoutScheme->GetRedoutTime(out_idx); // time passes since application of the RF pulse float tRf = m_Parameters->m_SignalGen.m_tEcho+t; // calculate eddy current decay factor // (TODO: vielleicht umbauen dass hier die zeit vom letzten diffusionsgradienten an genommen wird. doku dann auch entsprechend anpassen.) float eddyDecay = 0; if ( m_Parameters->m_Misc.m_DoAddEddyCurrents && m_Parameters->m_SignalGen.m_EddyStrength>0) eddyDecay = std::exp(-tRead/m_Parameters->m_SignalGen.m_Tau ); // calcualte signal relaxation factors std::vector< float > relaxFactor; if ( m_Parameters->m_SignalGen.m_DoSimulateRelaxation) for (unsigned int i=0; im_SignalGen.m_tInhom)); } // add ghosting by adding gradient delay induced offset if (m_Parameters->m_Misc.m_DoAddGhosts) { if (out_idx[1]%2 == 1) kx -= m_Parameters->m_SignalGen.m_KspaceLineOffset; else kx += m_Parameters->m_SignalGen.m_KspaceLineOffset; } // pull stuff out of inner loop t /= 1000; kx /= xMax; ky /= yMaxFov; // calculate signal s at k-space position (kx, ky) vcl_complex s(0,0); InputIteratorType it(m_CompartmentImages[0], m_CompartmentImages[0]->GetLargestPossibleRegion() ); while( !it.IsAtEnd() ) { typename InputImageType::IndexType input_idx = it.GetIndex(); // shift x,y for DFT: (0 -- N) --> (-N/2 -- N/2) float x = input_idx[0]; float y = input_idx[1]; if (static_cast(xMax)%2==1) x -= (xMax-1)/2; else x -= xMax/2; if (static_cast(yMax)%2==1) y -= (yMax-1)/2; else y -= yMax/2; // sum compartment signals and simulate relaxation ScalarType f_real = 0; for (unsigned int i=0; im_SignalGen.m_DoSimulateRelaxation) f_real += m_CompartmentImages[i]->GetPixel(input_idx) * relaxFactor[i] * m_Parameters->m_SignalGen.m_SignalScale; else f_real += m_CompartmentImages[i]->GetPixel(input_idx) * m_Parameters->m_SignalGen.m_SignalScale; // vector from image center to current position (in meter) // only necessary for eddy currents and non-constant coil sensitivity VectorType pos; if ((m_Parameters->m_Misc.m_DoAddEddyCurrents && m_Parameters->m_SignalGen.m_EddyStrength>0 && !m_IsBaseline) || m_Parameters->m_SignalGen.m_CoilSensitivityProfile!=SignalGenerationParameters::COIL_CONSTANT) { pos[0] = x; pos[1] = y; pos[2] = m_Z; pos = m_Transform*pos; } if (m_Parameters->m_SignalGen.m_CoilSensitivityProfile!=SignalGenerationParameters::COIL_CONSTANT) f_real *= CoilSensitivity(pos); // simulate eddy currents and other distortions float omega = 0; // frequency offset if ( m_Parameters->m_Misc.m_DoAddEddyCurrents && m_Parameters->m_SignalGen.m_EddyStrength>0 && !m_IsBaseline) omega += (m_DiffusionGradientDirection[0]*pos[0]+m_DiffusionGradientDirection[1]*pos[1]+m_DiffusionGradientDirection[2]*pos[2]) * eddyDecay; // simulate distortions if (m_Parameters->m_Misc.m_DoAddDistortions) { if (m_MovedFmap.IsNotNull()) // if we have headmotion, use moved map omega += m_MovedFmap->GetPixel(input_idx); else if (m_Parameters->m_SignalGen.m_FrequencyMap.IsNotNull()) { itk::Image::IndexType index; index[0] = input_idx[0]; index[1] = input_idx[1]; index[2] = m_Zidx; omega += m_Parameters->m_SignalGen.m_FrequencyMap->GetPixel(index); } } // if signal comes from outside FOV, mirror it back (wrap-around artifact - aliasing if (m_Parameters->m_Misc.m_DoAddAliasing) { if (y<-yMaxFov_half) y += yMaxFov; else if (y>=yMaxFov_half) y -= yMaxFov; } // actual DFT term vcl_complex f(f_real, 0); s += f * std::exp( std::complex(0, itk::Math::twopi * (kx*x + ky*y + omega*t )) ); ++it; } 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; s += outputImage->GetPixel(kIdx); // add precalculated noise outputImage->SetPixel(kIdx, s); m_KSpaceImage->SetPixel(kIdx, sqrt(s.imag()*s.imag()+s.real()*s.real()) ); ++oit; } } template< class ScalarType > void KspaceImageFilter< ScalarType > ::AfterThreadedGenerateData() { delete m_ReadoutScheme; typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); int kxMax = outputImage->GetLargestPossibleRegion().GetSize(0); // k-space size in x-direction int 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]; // reverse phase if (!m_Parameters->m_SignalGen.m_ReversePhase) kIdx[1] = static_cast(kyMax-1-kIdx[1]); if (kIdx[1]>kyMax*m_Parameters->m_SignalGen.m_PartialFourier) { // reverse readout direction if (oit.GetIndex()[1]%2 == 1) kIdx[0] = kxMax-kIdx[0]-1; // calculate symmetric index itk::Index< 2 > kIdx2; kIdx2[0] = (kxMax-kIdx[0]-kxMax%2)%kxMax; kIdx2[1] = (kyMax-kIdx[1]-kyMax%2)%kyMax; // use complex conjugate of symmetric index value at current index vcl_complex s = outputImage->GetPixel(kIdx2); s = vcl_complex(s.real(), -s.imag()); outputImage->SetPixel(kIdx, s); m_KSpaceImage->SetPixel(kIdx, 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()%kxMax; spikeIdx[1] = m_RandGen->GetIntegerVariate()%kyMax; outputImage->SetPixel(spikeIdx, m_Spike); m_SpikeLog += "[" + boost::lexical_cast(spikeIdx[0]) + "," + boost::lexical_cast(spikeIdx[1]) + "," + boost::lexical_cast(m_Zidx) + "] Magnitude: " + boost::lexical_cast(m_Spike.real()) + "+" + boost::lexical_cast(m_Spike.imag()) + "i\n"; } } } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp index c81c025fa3..86a0959de0 100755 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp @@ -1,1790 +1,1790 @@ /*=================================================================== 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 #include #include #include #include #include #include #include #include #include namespace itk { template< class PixelType > TractsToDWIImageFilter< PixelType >::TractsToDWIImageFilter() : m_StatusText("") , m_UseConstantRandSeed(false) , m_RandGen(itk::Statistics::MersenneTwisterRandomVariateGenerator::New()) { m_DoubleInterpolator = itk::LinearInterpolateImageFunction< ItkDoubleImgType, float >::New(); m_NullDir.Fill(0); } template< class PixelType > TractsToDWIImageFilter< PixelType >::~TractsToDWIImageFilter() { } template< class PixelType > TractsToDWIImageFilter< PixelType >::DoubleDwiType::Pointer TractsToDWIImageFilter< PixelType >:: SimulateKspaceAcquisition( std::vector< DoubleDwiType::Pointer >& compartment_images ) { unsigned int numFiberCompartments = m_Parameters.m_FiberModelList.size(); // create slice object ImageRegion<2> sliceRegion; sliceRegion.SetSize(0, m_WorkingImageRegion.GetSize()[0]); sliceRegion.SetSize(1, m_WorkingImageRegion.GetSize()[1]); Vector< double, 2 > sliceSpacing; sliceSpacing[0] = m_WorkingSpacing[0]; sliceSpacing[1] = m_WorkingSpacing[1]; DoubleDwiType::PixelType nullPix; nullPix.SetSize(compartment_images.at(0)->GetVectorLength()); nullPix.Fill(0.0); auto magnitudeDwiImage = DoubleDwiType::New(); magnitudeDwiImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); magnitudeDwiImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin ); magnitudeDwiImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); magnitudeDwiImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); magnitudeDwiImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); magnitudeDwiImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); magnitudeDwiImage->SetVectorLength( compartment_images.at(0)->GetVectorLength() ); magnitudeDwiImage->Allocate(); magnitudeDwiImage->FillBuffer(nullPix); m_PhaseImage = DoubleDwiType::New(); m_PhaseImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); m_PhaseImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin ); m_PhaseImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); m_PhaseImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_PhaseImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_PhaseImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_PhaseImage->SetVectorLength( compartment_images.at(0)->GetVectorLength() ); m_PhaseImage->Allocate(); m_PhaseImage->FillBuffer(nullPix); m_KspaceImage = DoubleDwiType::New(); m_KspaceImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); m_KspaceImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin ); m_KspaceImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); m_KspaceImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_KspaceImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_KspaceImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_KspaceImage->SetVectorLength( m_Parameters.m_SignalGen.m_NumberOfCoils ); m_KspaceImage->Allocate(); m_KspaceImage->FillBuffer(nullPix); // calculate coil positions double a = m_Parameters.m_SignalGen.m_ImageRegion.GetSize(0)*m_Parameters.m_SignalGen.m_ImageSpacing[0]; double b = m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1)*m_Parameters.m_SignalGen.m_ImageSpacing[1]; double c = m_Parameters.m_SignalGen.m_ImageRegion.GetSize(2)*m_Parameters.m_SignalGen.m_ImageSpacing[2]; double diagonal = sqrt(a*a+b*b)/1000; // image diagonal in m m_CoilPointset = mitk::PointSet::New(); std::vector< itk::Vector > coilPositions; itk::Vector pos; pos.Fill(0.0); pos[1] = -diagonal/2; itk::Vector center; center[0] = a/2-m_Parameters.m_SignalGen.m_ImageSpacing[0]/2; center[1] = b/2-m_Parameters.m_SignalGen.m_ImageSpacing[2]/2; center[2] = c/2-m_Parameters.m_SignalGen.m_ImageSpacing[1]/2; for (unsigned int c=0; cInsertPoint(c, pos*1000 + m_Parameters.m_SignalGen.m_ImageOrigin.GetVectorFromOrigin() + center ); double rz = 360.0/m_Parameters.m_SignalGen.m_NumberOfCoils * itk::Math::pi/180; vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); rotZ[0][0] = cos(rz); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(rz); rotZ[1][0] = -rotZ[0][1]; pos.SetVnlVector(rotZ*pos.GetVnlVector()); } auto num_slices = compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2); auto num_gradient_volumes = static_cast(compartment_images.at(0)->GetVectorLength()); auto max_threads = omp_get_max_threads(); int out_threads = Math::ceil(std::sqrt(max_threads)); int in_threads = Math::floor(std::sqrt(max_threads)); if (out_threads > num_gradient_volumes) { out_threads = num_gradient_volumes; in_threads = Math::floor(static_cast(max_threads/out_threads)); } PrintToLog("Parallel volumes: " + boost::lexical_cast(out_threads), false, true, true); PrintToLog("Threads per slice: " + boost::lexical_cast(in_threads), false, true, true); std::list< std::tuple > spikes; if (m_Parameters.m_Misc.m_DoAddSpikes) for (unsigned int i=0; i( m_RandGen->GetIntegerVariate()%num_gradient_volumes, m_RandGen->GetIntegerVariate()%num_slices, m_RandGen->GetIntegerVariate()%m_Parameters.m_SignalGen.m_NumberOfCoils); spikes.push_back(spike); } PrintToLog("0% 10 20 30 40 50 60 70 80 90 100%", false, true, false); PrintToLog("|----|----|----|----|----|----|----|----|----|----|\n*", false, false, false); unsigned long lastTick = 0; boost::progress_display disp(static_cast(num_gradient_volumes)*compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2)); #pragma omp parallel for num_threads(out_threads) for (int g=0; gGetAbortGenerateData()) continue; std::list< std::tuple > spikeSlice; #pragma omp critical { for (auto spike : spikes) if (std::get<0>(spike) == static_cast(g)) spikeSlice.push_back(std::tuple(std::get<1>(spike), std::get<2>(spike))); } for (unsigned int z=0; z compartment_slices; std::vector< float > t2Vector; std::vector< float > 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++) { Float2DImageType::IndexType index2D; index2D[0]=x; index2D[1]=y; DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z; slice->SetPixel(index2D, compartment_images.at(i)->GetPixel(index3D)[g]); } compartment_slices.push_back(slice); t2Vector.push_back(signalModel->GetT2()); t1Vector.push_back(signalModel->GetT1()); } if (this->GetAbortGenerateData()) continue; for (unsigned int c=0; c(ss) == z && std::get<1>(ss) == c) ++numSpikes; // create k-sapce (inverse fourier transform slices) auto idft = itk::KspaceImageFilter< Float2DImageType::PixelType >::New(); idft->SetCompartmentImages(compartment_slices); idft->SetT2(t2Vector); idft->SetT1(t1Vector); if (m_UseConstantRandSeed) { int linear_seed = g + num_gradient_volumes*z + num_gradient_volumes*compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2)*c; idft->SetRandSeed(linear_seed); } idft->SetParameters(&m_Parameters); idft->SetZ((float)z-(float)( compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2) -compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2)%2 ) / 2.0); idft->SetZidx(z); idft->SetCoilPosition(coilPositions.at(c)); idft->SetFiberBundle(m_FiberBundle); idft->SetTranslation(m_Translations.at(g)); idft->SetRotationMatrix(m_RotationsInv.at(g)); idft->SetDiffusionGradientDirection(m_Parameters.m_SignalGen.GetGradientDirection(g)); idft->SetSpikesPerSlice(numSpikes); idft->SetNumberOfThreads(in_threads); idft->Update(); #pragma omp critical if (numSpikes>0) { m_SpikeLog += "Volume " + boost::lexical_cast(g) + " Coil " + boost::lexical_cast(c) + "\n"; m_SpikeLog += idft->GetSpikeLog(); } Complex2DImageType::Pointer fSlice; fSlice = idft->GetOutput(); // fourier transform slice Complex2DImageType::Pointer newSlice; auto dft = itk::DftImageFilter< Float2DImageType::PixelType >::New(); dft->SetInput(fSlice); dft->SetParameters(m_Parameters); dft->SetNumberOfThreads(in_threads); 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; Complex2DImageType::IndexType index2D; index2D[0]=x; index2D[1]=y; Complex2DImageType::PixelType cPix = newSlice->GetPixel(index2D); double magn = sqrt(cPix.real()*cPix.real()+cPix.imag()*cPix.imag()); double phase = 0; if (cPix.real()!=0) phase = atan( cPix.imag()/cPix.real() ); DoubleDwiType::PixelType real_pix = m_OutputImagesReal.at(c)->GetPixel(index3D); real_pix[g] = cPix.real(); m_OutputImagesReal.at(c)->SetPixel(index3D, real_pix); DoubleDwiType::PixelType imag_pix = m_OutputImagesImag.at(c)->GetPixel(index3D); imag_pix[g] = cPix.imag(); m_OutputImagesImag.at(c)->SetPixel(index3D, imag_pix); DoubleDwiType::PixelType dwiPix = magnitudeDwiImage->GetPixel(index3D); DoubleDwiType::PixelType phasePix = m_PhaseImage->GetPixel(index3D); if (m_Parameters.m_SignalGen.m_NumberOfCoils>1) { dwiPix[g] += magn*magn; phasePix[g] += phase*phase; } else { dwiPix[g] = magn; phasePix[g] = phase; } //#pragma omp critical { magnitudeDwiImage->SetPixel(index3D, dwiPix); m_PhaseImage->SetPixel(index3D, phasePix); // k-space image if (g==0) { DoubleDwiType::PixelType kspacePix = m_KspaceImage->GetPixel(index3D); kspacePix[c] = idft->GetKSpaceImage()->GetPixel(index2D); m_KspaceImage->SetPixel(index3D, kspacePix); } } } } if (m_Parameters.m_SignalGen.m_NumberOfCoils>1) { for (int y=0; y(magnitudeDwiImage->GetLargestPossibleRegion().GetSize(1)); y++) for (int x=0; x(magnitudeDwiImage->GetLargestPossibleRegion().GetSize(0)); x++) { DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z; DoubleDwiType::PixelType magPix = magnitudeDwiImage->GetPixel(index3D); magPix[g] = sqrt(magPix[g]/m_Parameters.m_SignalGen.m_NumberOfCoils); DoubleDwiType::PixelType phasePix = m_PhaseImage->GetPixel(index3D); phasePix[g] = sqrt(phasePix[g]/m_Parameters.m_SignalGen.m_NumberOfCoils); //#pragma omp critical { magnitudeDwiImage->SetPixel(index3D, magPix); m_PhaseImage->SetPixel(index3D, phasePix); } } } ++disp; unsigned long newTick = 50*disp.count()/disp.expected_count(); for (unsigned long tick = 0; tick<(newTick-lastTick); tick++) PrintToLog("*", false, false, false); lastTick = newTick; } } PrintToLog("\n", false); return magnitudeDwiImage; } template< class PixelType > TractsToDWIImageFilter< PixelType >::ItkDoubleImgType::Pointer TractsToDWIImageFilter< PixelType >:: NormalizeInsideMask(ItkDoubleImgType::Pointer image) { double max = itk::NumericTraits< double >::min(); double min = itk::NumericTraits< double >::max(); itk::ImageRegionIterator< ItkDoubleImgType > it(image, image->GetLargestPossibleRegion()); while(!it.IsAtEnd()) { if (m_Parameters.m_SignalGen.m_MaskImage.IsNotNull() && m_Parameters.m_SignalGen.m_MaskImage->GetPixel(it.GetIndex())<=0) { it.Set(0.0); ++it; continue; } if (it.Get()>max) max = it.Get(); if (it.Get()::New(); scaler->SetInput(image); scaler->SetShift(-min); scaler->SetScale(1.0/(max-min)); scaler->Update(); return scaler->GetOutput(); } template< class PixelType > void TractsToDWIImageFilter< PixelType >::CheckVolumeFractionImages() { m_UseRelativeNonFiberVolumeFractions = false; // check for fiber volume fraction maps unsigned int fibVolImages = 0; for (std::size_t i=0; iGetVolumeFractionImage().IsNotNull()) { PrintToLog("Using volume fraction map for fiber compartment " + boost::lexical_cast(i+1), false); fibVolImages++; } } // check for non-fiber volume fraction maps unsigned int nonfibVolImages = 0; for (std::size_t i=0; iGetVolumeFractionImage().IsNotNull()) { PrintToLog("Using volume fraction map for non-fiber compartment " + boost::lexical_cast(i+1), false); nonfibVolImages++; } } // not all fiber compartments are using volume fraction maps // --> non-fiber volume fractions are assumed to be relative to the // non-fiber volume and not absolute voxel-volume fractions. // this means if two non-fiber compartments are used but only one of them // has an associated volume fraction map, the repesctive other volume fraction map // can be determined as inverse (1-val) of the present volume fraction map- if ( fibVolImages::New(); inverter->SetMaximum(1.0); if ( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage().IsNull() && m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage().IsNotNull() ) { // m_Parameters.m_NonFiberModelList[1]->SetVolumeFractionImage( // NormalizeInsideMask( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage() ) ); inverter->SetInput( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage() ); inverter->Update(); m_Parameters.m_NonFiberModelList[0]->SetVolumeFractionImage(inverter->GetOutput()); } else if ( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage().IsNull() && m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage().IsNotNull() ) { // m_Parameters.m_NonFiberModelList[0]->SetVolumeFractionImage( // NormalizeInsideMask( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage() ) ); inverter->SetInput( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage() ); inverter->Update(); m_Parameters.m_NonFiberModelList[1]->SetVolumeFractionImage(inverter->GetOutput()); } else { itkExceptionMacro("Something went wrong in automatically calculating the missing non-fiber volume fraction image!" " Did you use two non fiber compartments but only one volume fraction image?" " Then it should work and this error is really strange."); } m_UseRelativeNonFiberVolumeFractions = true; nonfibVolImages++; } // Up to two fiber compartments are allowed without volume fraction maps since the volume fractions can then be determined automatically if (m_Parameters.m_FiberModelList.size()>2 && fibVolImages!=m_Parameters.m_FiberModelList.size()) itkExceptionMacro("More than two fiber compartment selected but no corresponding volume fraction maps set!"); // One non-fiber compartment is allowed without volume fraction map since the volume fraction can then be determined automatically if (m_Parameters.m_NonFiberModelList.size()>1 && nonfibVolImages!=m_Parameters.m_NonFiberModelList.size()) itkExceptionMacro("More than one non-fiber compartment selected but no volume fraction maps set!"); if (fibVolImages0) { PrintToLog("Not all fiber compartments are using an associated volume fraction image.\n" "Assuming non-fiber volume fraction images to contain values relative to the" " remaining non-fiber volume, not absolute values.", false); m_UseRelativeNonFiberVolumeFractions = true; // mitk::LocaleSwitch localeSwitch("C"); // itk::ImageFileWriter::Pointer wr = itk::ImageFileWriter::New(); // wr->SetInput(m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage()); // wr->SetFileName("/local/volumefraction.nrrd"); // wr->Update(); } // initialize the images that store the output volume fraction of each compartment m_VolumeFractions.clear(); for (std::size_t i=0; 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() { m_Rotations.clear(); m_Translations.clear(); m_MotionLog = ""; m_SpikeLog = ""; // initialize output dwi image m_Parameters.m_SignalGen.m_CroppedRegion = m_Parameters.m_SignalGen.m_ImageRegion; if (m_Parameters.m_Misc.m_DoAddAliasing) m_Parameters.m_SignalGen.m_CroppedRegion.SetSize( 1, m_Parameters.m_SignalGen.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_Parameters.m_SignalGen.m_CroppedRegion.GetSize(1))*m_Parameters.m_SignalGen.m_ImageSpacing[1]/2; m_OutputImage = OutputImageType::New(); m_OutputImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); m_OutputImage->SetOrigin( shiftedOrigin ); m_OutputImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); m_OutputImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_OutputImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_OutputImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_OutputImage->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() ); m_OutputImage->Allocate(); typename OutputImageType::PixelType temp; temp.SetSize(m_Parameters.m_SignalGen.GetNumVolumes()); temp.Fill(0.0); m_OutputImage->FillBuffer(temp); PrintToLog("Output image spacing: [" + boost::lexical_cast(m_Parameters.m_SignalGen.m_ImageSpacing[0]) + "," + boost::lexical_cast(m_Parameters.m_SignalGen.m_ImageSpacing[1]) + "," + boost::lexical_cast(m_Parameters.m_SignalGen.m_ImageSpacing[2]) + "]", false); PrintToLog("Output image size: [" + boost::lexical_cast(m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(0)) + "," + boost::lexical_cast(m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(1)) + "," + boost::lexical_cast(m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(2)) + "]", false); // images containing real and imaginary part of the dMRI signal for each coil m_OutputImagesReal.clear(); m_OutputImagesImag.clear(); for (unsigned int i=0; iSetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); outputImageReal->SetOrigin( shiftedOrigin ); outputImageReal->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); outputImageReal->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); outputImageReal->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); outputImageReal->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); outputImageReal->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() ); outputImageReal->Allocate(); outputImageReal->FillBuffer(temp); m_OutputImagesReal.push_back(outputImageReal); typename DoubleDwiType::Pointer outputImageImag = DoubleDwiType::New(); outputImageImag->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); outputImageImag->SetOrigin( shiftedOrigin ); outputImageImag->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); outputImageImag->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); outputImageImag->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); outputImageImag->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); outputImageImag->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() ); outputImageImag->Allocate(); outputImageImag->FillBuffer(temp); m_OutputImagesImag.push_back(outputImageImag); } // Apply in-plane upsampling for Gibbs ringing artifact double upsampling = 1; if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging && m_Parameters.m_SignalGen.m_ZeroRinging==0) 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]; PrintToLog("Working image spacing: [" + boost::lexical_cast(m_WorkingSpacing[0]) + "," + boost::lexical_cast(m_WorkingSpacing[1]) + "," + boost::lexical_cast(m_WorkingSpacing[2]) + "]", false); PrintToLog("Working image size: [" + boost::lexical_cast(m_WorkingImageRegion.GetSize(0)) + "," + boost::lexical_cast(m_WorkingImageRegion.GetSize(1)) + "," + boost::lexical_cast(m_WorkingImageRegion.GetSize(2)) + "]", false); // 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); } if (m_FiberBundle.IsNull() && m_InputImage.IsNotNull()) { m_CompartmentImages.clear(); m_Parameters.m_SignalGen.m_DoAddMotion = false; m_Parameters.m_SignalGen.m_DoSimulateRelaxation = false; PrintToLog("Simulating acquisition for input diffusion-weighted image.", false); auto caster = itk::CastImageFilter< OutputImageType, DoubleDwiType >::New(); caster->SetInput(m_InputImage); caster->Update(); if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging && m_Parameters.m_SignalGen.m_ZeroRinging==0) { PrintToLog("Upsampling input diffusion-weighted image for Gibbs ringing simulation.", false); auto resampler = itk::ResampleDwiImageFilter< double >::New(); resampler->SetInput(caster->GetOutput()); itk::Vector< double, 3 > samplingFactor; samplingFactor[0] = upsampling; samplingFactor[1] = upsampling; samplingFactor[2] = 1; resampler->SetSamplingFactor(samplingFactor); resampler->SetInterpolation(itk::ResampleDwiImageFilter< double >::Interpolate_WindowedSinc); resampler->Update(); m_CompartmentImages.push_back(resampler->GetOutput()); } else m_CompartmentImages.push_back(caster->GetOutput()); VectorType translation; translation.Fill(0.0); MatrixType rotation; rotation.SetIdentity(); for (unsigned int g=0; gGetLargestPossibleRegion()!=m_WorkingImageRegion) { PrintToLog("Resampling tissue mask", false); // rescale mask image (otherwise there are problems with the resampling) auto rescaler = itk::RescaleIntensityImageFilter::New(); rescaler->SetInput(0,m_Parameters.m_SignalGen.m_MaskImage); rescaler->SetOutputMaximum(100); rescaler->SetOutputMinimum(0); rescaler->Update(); // resample mask image auto resampler = itk::ResampleImageFilter::New(); resampler->SetInput(rescaler->GetOutput()); resampler->SetSize(m_WorkingImageRegion.GetSize()); resampler->SetOutputSpacing(m_WorkingSpacing); resampler->SetOutputOrigin(m_WorkingOrigin); resampler->SetOutputDirection(m_Parameters.m_SignalGen.m_ImageDirection); resampler->SetOutputStartIndex ( m_WorkingImageRegion.GetIndex() ); auto nn_interpolator = itk::NearestNeighborInterpolateImageFunction::New(); resampler->SetInterpolator(nn_interpolator); resampler->Update(); m_Parameters.m_SignalGen.m_MaskImage = resampler->GetOutput(); } // resample frequency map if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull() && m_Parameters.m_SignalGen.m_FrequencyMap->GetLargestPossibleRegion()!=m_WorkingImageRegion) { PrintToLog("Resampling frequency map", false); auto resampler = itk::ResampleImageFilter::New(); resampler->SetInput(m_Parameters.m_SignalGen.m_FrequencyMap); resampler->SetSize(m_WorkingImageRegion.GetSize()); resampler->SetOutputSpacing(m_WorkingSpacing); resampler->SetOutputOrigin(m_WorkingOrigin); resampler->SetOutputDirection(m_Parameters.m_SignalGen.m_ImageDirection); resampler->SetOutputStartIndex ( m_WorkingImageRegion.GetIndex() ); auto nn_interpolator = itk::NearestNeighborInterpolateImageFunction::New(); resampler->SetInterpolator(nn_interpolator); resampler->Update(); m_Parameters.m_SignalGen.m_FrequencyMap = resampler->GetOutput(); } m_MaskImageSet = true; if (m_Parameters.m_SignalGen.m_MaskImage.IsNull()) { // no input tissue mask is set -> create default PrintToLog("No tissue mask set", false); m_Parameters.m_SignalGen.m_MaskImage = ItkUcharImgType::New(); m_Parameters.m_SignalGen.m_MaskImage->SetSpacing( m_WorkingSpacing ); m_Parameters.m_SignalGen.m_MaskImage->SetOrigin( m_WorkingOrigin ); m_Parameters.m_SignalGen.m_MaskImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); m_Parameters.m_SignalGen.m_MaskImage->SetLargestPossibleRegion( m_WorkingImageRegion ); m_Parameters.m_SignalGen.m_MaskImage->SetBufferedRegion( m_WorkingImageRegion ); m_Parameters.m_SignalGen.m_MaskImage->SetRequestedRegion( m_WorkingImageRegion ); m_Parameters.m_SignalGen.m_MaskImage->Allocate(); m_Parameters.m_SignalGen.m_MaskImage->FillBuffer(100); m_MaskImageSet = false; } else { PrintToLog("Using tissue mask", false); } if (m_Parameters.m_SignalGen.m_DoAddMotion) { if (m_Parameters.m_SignalGen.m_DoRandomizeMotion) { PrintToLog("Random motion artifacts:", false); PrintToLog("Maximum rotation: +/-" + boost::lexical_cast(m_Parameters.m_SignalGen.m_Rotation) + "°", false); PrintToLog("Maximum translation: +/-" + boost::lexical_cast(m_Parameters.m_SignalGen.m_Translation) + "mm", false); } else { PrintToLog("Linear motion artifacts:", false); PrintToLog("Maximum rotation: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_Rotation) + "°", false); PrintToLog("Maximum translation: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_Translation) + "mm", false); } } if ( m_Parameters.m_SignalGen.m_MotionVolumes.empty() ) { // no motion in first volume m_Parameters.m_SignalGen.m_MotionVolumes.push_back(false); // motion in all other volumes while ( m_Parameters.m_SignalGen.m_MotionVolumes.size() < m_Parameters.m_SignalGen.GetNumVolumes() ) { m_Parameters.m_SignalGen.m_MotionVolumes.push_back(true); } } // we need to know for every volume if there is motion. if this information is missing, then set corresponding fal to false while ( m_Parameters.m_SignalGen.m_MotionVolumes.size()::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; VectorType 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(); auto 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); auto nn_interpolator = itk::NearestNeighborInterpolateImageFunction::New(); upsampler->SetInterpolator(nn_interpolator); upsampler->Update(); m_UpsampledMaskImage = upsampler->GetOutput(); } template< class PixelType > void TractsToDWIImageFilter< PixelType >::InitializeFiberData() { m_mmRadius = m_Parameters.m_SignalGen.m_AxonRadius/1000; auto caster = itk::CastImageFilter< itk::Image, itk::Image >::New(); caster->SetInput(m_TransformedMaskImage); caster->Update(); vtkSmartPointer weights = m_FiberBundle->GetFiberWeights(); float mean_weight = 0; for (int i=0; iGetSize(); i++) mean_weight += weights->GetValue(i); mean_weight /= weights->GetSize(); if (mean_weight>0.000001) for (int i=0; iGetSize(); i++) m_FiberBundle->SetFiberWeight(i, weights->GetValue(i)/mean_weight); else PrintToLog("\nWarning: streamlines have VERY low weights. Average weight: " + boost::lexical_cast(mean_weight) + ". Possible source of calculation errors.", false, true, true); auto density_calculator = itk::TractDensityImageFilter< itk::Image >::New(); density_calculator->SetFiberBundle(m_FiberBundle); density_calculator->SetInputImage(caster->GetOutput()); density_calculator->SetBinaryOutput(false); density_calculator->SetUseImageGeometry(true); density_calculator->SetOutputAbsoluteValues(true); density_calculator->Update(); double max_density = density_calculator->GetMaxDensity(); double voxel_volume = m_WorkingSpacing[0]*m_WorkingSpacing[1]*m_WorkingSpacing[2]; if (m_mmRadius>0) { std::stringstream stream; stream << std::fixed << setprecision(2) << itk::Math::pi*m_mmRadius*m_mmRadius*max_density; std::string s = stream.str(); PrintToLog("\nMax. fiber volume: " + s + "mm².", false, true, true); { double full_radius = 1000*std::sqrt(voxel_volume/(max_density*itk::Math::pi)); std::stringstream stream; stream << std::fixed << setprecision(2) << full_radius; std::string s = stream.str(); PrintToLog("\nA full fiber voxel corresponds to a fiber radius of ~" + s + "µm, given the current fiber configuration.", false, true, true); } } else { m_mmRadius = std::sqrt(voxel_volume/(max_density*itk::Math::pi)); std::stringstream stream; stream << std::fixed << setprecision(2) << m_mmRadius*1000; std::string s = stream.str(); PrintToLog("\nSetting fiber radius to " + s + "µm to obtain full voxel.", false, true, true); } // a second fiber bundle is needed to store the transformed version of the m_FiberBundleWorkingCopy m_FiberBundleTransformed = m_FiberBundle->GetDeepCopy(); } template< class PixelType > bool TractsToDWIImageFilter< PixelType >::PrepareLogFile() { if(m_Logfile.is_open()) m_Logfile.close(); std::string filePath; std::string fileName; // Get directory name: if (m_Parameters.m_Misc.m_OutputPath.size() > 0) { filePath = m_Parameters.m_Misc.m_OutputPath; if( *(--(filePath.cend())) != '/') { filePath.push_back('/'); } } else { filePath = mitk::IOUtil::GetTempPath() + '/'; } // check if directory exists, else use /tmp/: if( itksys::SystemTools::FileIsDirectory( filePath ) ) { while( *(--(filePath.cend())) == '/') { filePath.pop_back(); } filePath = filePath + '/'; } else { filePath = mitk::IOUtil::GetTempPath() + '/'; } // Get file name: if( ! m_Parameters.m_Misc.m_ResultNode->GetName().empty() ) { fileName = m_Parameters.m_Misc.m_ResultNode->GetName(); } else { fileName = ""; } if( ! m_Parameters.m_Misc.m_OutputPrefix.empty() ) { fileName = m_Parameters.m_Misc.m_OutputPrefix + fileName; } else { fileName = "fiberfox"; } // check if file already exists and DO NOT overwrite existing files: std::string NameTest = fileName; int c = 0; while( itksys::SystemTools::FileExists( filePath + '/' + fileName + ".log" ) && c <= std::numeric_limits::max() ) { fileName = NameTest + "_" + boost::lexical_cast(c); ++c; } try { m_Logfile.open( ( filePath + '/' + fileName + ".log" ).c_str() ); } catch (const std::ios_base::failure &fail) { MITK_ERROR << "itkTractsToDWIImageFilter.cpp: Exception " << fail.what() << " while trying to open file" << filePath << '/' << fileName << ".log"; return false; } if ( m_Logfile.is_open() ) { PrintToLog( "Logfile: " + filePath + '/' + fileName + ".log", false ); return true; } else { m_StatusText += "Logfile could not be opened!\n"; MITK_ERROR << "itkTractsToDWIImageFilter.cpp: Logfile could not be opened!"; return false; } } template< class PixelType > void TractsToDWIImageFilter< PixelType >::GenerateData() { PrintToLog("\n**********************************************", false); // prepare logfile if ( ! PrepareLogFile() ) { this->SetAbortGenerateData( true ); return; } PrintToLog("Starting Fiberfox dMRI simulation"); m_TimeProbe.Start(); // check input data if (m_FiberBundle.IsNull() && m_InputImage.IsNull()) itkExceptionMacro("Input fiber bundle and input diffusion-weighted image is nullptr!"); if (m_Parameters.m_FiberModelList.empty() && m_InputImage.IsNull()) itkExceptionMacro("No diffusion model for fiber compartments defined and input diffusion-weighted" " image is nullptr! At least one fiber compartment is necessary to simulate diffusion."); if (m_Parameters.m_NonFiberModelList.empty() && m_InputImage.IsNull()) itkExceptionMacro("No diffusion model for non-fiber compartments defined and input diffusion-weighted" " image is nullptr! At least one non-fiber compartment is necessary to simulate diffusion."); if (m_Parameters.m_SignalGen.m_DoDisablePartialVolume) // no partial volume? remove all but first fiber compartment while (m_Parameters.m_FiberModelList.size()>1) m_Parameters.m_FiberModelList.pop_back(); if (!m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition) // No upsampling of input image needed if no k-space simulation is performed m_Parameters.m_SignalGen.m_DoAddGibbsRinging = false; if (m_UseConstantRandSeed) // always generate the same random numbers? m_RandGen->SetSeed(0); else m_RandGen->SetSeed(); InitializeData(); if ( m_FiberBundle.IsNotNull() ) // if no fiber bundle is found, we directly proceed to the k-space acquisition simulation { CheckVolumeFractionImages(); InitializeFiberData(); int numFiberCompartments = m_Parameters.m_FiberModelList.size(); int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size(); double maxVolume = 0; unsigned long lastTick = 0; int signalModelSeed = m_RandGen->GetIntegerVariate(); PrintToLog("\n", false, false); PrintToLog("Generating " + boost::lexical_cast(numFiberCompartments+numNonFiberCompartments) + "-compartment diffusion-weighted signal."); std::vector< int > bVals = m_Parameters.m_SignalGen.GetBvalues(); PrintToLog("b-values: ", false, false, true); for (auto v : bVals) PrintToLog(boost::lexical_cast(v) + " ", false, false, true); PrintToLog("\nVolumes: " + boost::lexical_cast(m_Parameters.m_SignalGen.GetNumVolumes()), false, true, true); PrintToLog("\n", false, false, true); PrintToLog("\n", false, false, true); unsigned int image_size_x = m_WorkingImageRegion.GetSize(0); unsigned int region_size_y = m_WorkingImageRegion.GetSize(1); unsigned int num_gradients = m_Parameters.m_SignalGen.GetNumVolumes(); int numFibers = m_FiberBundle->GetNumFibers(); boost::progress_display disp(numFibers*num_gradients); if (m_FiberBundle->GetMeanFiberLength()<5.0) omp_set_num_threads(2); PrintToLog("0% 10 20 30 40 50 60 70 80 90 100%", false, true, false); PrintToLog("|----|----|----|----|----|----|----|----|----|----|\n*", false, false, false); for (unsigned int g=0; gSetSeed(signalModelSeed); for (std::size_t i=0; iSetSeed(signalModelSeed); // storing voxel-wise intra-axonal volume in mm³ auto intraAxonalVolumeImage = ItkDoubleImgType::New(); intraAxonalVolumeImage->SetSpacing( m_WorkingSpacing ); intraAxonalVolumeImage->SetOrigin( m_WorkingOrigin ); intraAxonalVolumeImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); intraAxonalVolumeImage->SetLargestPossibleRegion( m_WorkingImageRegion ); intraAxonalVolumeImage->SetBufferedRegion( m_WorkingImageRegion ); intraAxonalVolumeImage->SetRequestedRegion( m_WorkingImageRegion ); intraAxonalVolumeImage->Allocate(); intraAxonalVolumeImage->FillBuffer(0); maxVolume = 0; double* intraAxBuffer = intraAxonalVolumeImage->GetBufferPointer(); if (this->GetAbortGenerateData()) continue; vtkPolyData* fiberPolyData = m_FiberBundleTransformed->GetFiberPolyData(); // generate fiber signal (if there are any fiber models present) if (!m_Parameters.m_FiberModelList.empty()) { std::vector< double* > buffers; for (unsigned int i=0; iGetBufferPointer()); #pragma omp parallel for for( int i=0; iGetAbortGenerateData()) continue; float fiberWeight = m_FiberBundleTransformed->GetFiberWeight(i); int numPoints = -1; std::vector< itk::Vector > points_copy; #pragma omp critical { vtkCell* cell = fiberPolyData->GetCell(i); numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j))); } if (numPoints<2) continue; double seg_volume = fiberWeight*itk::Math::pi*m_mmRadius*m_mmRadius; for( int j=0; jGetAbortGenerateData()) { j=numPoints; continue; } itk::Vector v = points_copy.at(j); itk::Vector dir = points_copy.at(j+1)-v; if ( dir.GetSquaredNorm()<0.0001 || dir[0]!=dir[0] || dir[1]!=dir[1] || dir[2]!=dir[2] ) continue; dir.Normalize(); itk::Point startVertex = points_copy.at(j); itk::Index<3> startIndex; itk::ContinuousIndex startIndexCont; m_TransformedMaskImage->TransformPhysicalPointToIndex(startVertex, startIndex); m_TransformedMaskImage->TransformPhysicalPointToContinuousIndex(startVertex, startIndexCont); itk::Point endVertex = points_copy.at(j+1); itk::Index<3> endIndex; itk::ContinuousIndex endIndexCont; m_TransformedMaskImage->TransformPhysicalPointToIndex(endVertex, endIndex); m_TransformedMaskImage->TransformPhysicalPointToContinuousIndex(endVertex, endIndexCont); std::vector< std::pair< itk::Index<3>, double > > segments = mitk::imv::IntersectImage(m_WorkingSpacing, startIndex, endIndex, startIndexCont, endIndexCont); // generate signal for each fiber compartment for (int k=0; kSimulateMeasurement(g, dir)*seg_volume; for (std::pair< itk::Index<3>, double > seg : segments) { if (!m_TransformedMaskImage->GetLargestPossibleRegion().IsInside(seg.first) || m_TransformedMaskImage->GetPixel(seg.first)<=0) continue; double seg_signal = seg.second*signal_add; unsigned int linear_index = g + num_gradients*seg.first[0] + num_gradients*image_size_x*seg.first[1] + num_gradients*image_size_x*region_size_y*seg.first[2]; // update dMRI volume #pragma omp atomic buffers[k][linear_index] += seg_signal; // update fiber volume image if (k==0) { linear_index = seg.first[0] + image_size_x*seg.first[1] + image_size_x*region_size_y*seg.first[2]; #pragma omp atomic intraAxBuffer[linear_index] += seg.second*seg_volume; double vol = intraAxBuffer[linear_index]; if (vol>maxVolume) { maxVolume = vol; } } } } } #pragma omp critical { // progress report ++disp; unsigned long newTick = 50*disp.count()/disp.expected_count(); for (unsigned int tick = 0; tick<(newTick-lastTick); ++tick) PrintToLog("*", false, false, false); lastTick = newTick; } } } // axon radius not manually defined --> set fullest voxel (maxVolume) to full fiber voxel double density_correctiony_global = 1.0; if (m_Parameters.m_SignalGen.m_AxonRadius<0.0001) density_correctiony_global = m_VoxelVolume/maxVolume; // generate non-fiber signal ImageRegionIterator it3(m_TransformedMaskImage, m_TransformedMaskImage->GetLargestPossibleRegion()); while(!it3.IsAtEnd()) { if (it3.Get()>0) { DoubleDwiType::IndexType index = it3.GetIndex(); double iAxVolume = intraAxonalVolumeImage->GetPixel(index); // get non-transformed point (remove headmotion tranformation) // this point lives in the volume fraction image space itk::Point volume_fraction_point; if ( m_Parameters.m_SignalGen.m_DoAddMotion ) volume_fraction_point = GetMovedPoint(index, false); else m_TransformedMaskImage->TransformIndexToPhysicalPoint(index, volume_fraction_point); if (m_Parameters.m_SignalGen.m_DoDisablePartialVolume) { if (iAxVolume>0.0001) // scale fiber compartment to voxel { DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(index); pix[g] *= m_VoxelVolume/iAxVolume; m_CompartmentImages.at(0)->SetPixel(index, pix); if (g==0) m_VolumeFractions.at(0)->SetPixel(index, 1); } else { DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(index); pix[g] = 0; m_CompartmentImages.at(0)->SetPixel(index, pix); SimulateExtraAxonalSignal(index, volume_fraction_point, 0, g); } } else { // manually defined axon radius and voxel overflow --> rescale to voxel volume if ( m_Parameters.m_SignalGen.m_AxonRadius>=0.0001 && iAxVolume>m_VoxelVolume ) { for (int i=0; iGetPixel(index); pix[g] *= m_VoxelVolume/iAxVolume; m_CompartmentImages.at(i)->SetPixel(index, pix); } iAxVolume = m_VoxelVolume; } // if volume fraction image is set use it, otherwise use global scaling factor double density_correction_voxel = density_correctiony_global; if ( m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()!=nullptr && iAxVolume>0.0001 ) { m_DoubleInterpolator->SetInputImage(m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()); double volume_fraction = mitk::imv::GetImageValue(volume_fraction_point, true, m_DoubleInterpolator); if (volume_fraction<0) mitkThrow() << "Volume fraction image (index 1) contains negative values (intra-axonal compartment)!"; density_correction_voxel = m_VoxelVolume*volume_fraction/iAxVolume; // remove iAxVolume sclaing and scale to volume_fraction } else if (m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()!=nullptr) density_correction_voxel = 0.0; // adjust intra-axonal compartment volume by density correction factor DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(index); pix[g] *= density_correction_voxel; m_CompartmentImages.at(0)->SetPixel(index, pix); // normalize remaining fiber volume fractions (they are rescaled in SimulateExtraAxonalSignal) if (iAxVolume>0.0001) { for (int i=1; iGetPixel(index); pix[g] /= iAxVolume; m_CompartmentImages.at(i)->SetPixel(index, pix); } } else { for (int i=1; iGetPixel(index); pix[g] = 0; m_CompartmentImages.at(i)->SetPixel(index, pix); } } iAxVolume = density_correction_voxel*iAxVolume; // new intra-axonal volume = old intra-axonal volume * correction factor // simulate other compartments SimulateExtraAxonalSignal(index, volume_fraction_point, iAxVolume, g); } } ++it3; } } PrintToLog("\n", false); } if (this->GetAbortGenerateData()) { PrintToLog("\n", false, false); PrintToLog("Simulation aborted"); return; } DoubleDwiType::Pointer doubleOutImage; double signalScale = m_Parameters.m_SignalGen.m_SignalScale; if ( m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition ) // do k-space stuff { PrintToLog("\n", false, false); PrintToLog("Simulating k-space acquisition using " +boost::lexical_cast(m_Parameters.m_SignalGen.m_NumberOfCoils) +" coil(s)"); switch (m_Parameters.m_SignalGen.m_AcquisitionType) { case SignalGenerationParameters::SingleShotEpi: { PrintToLog("Acquisition type: single shot EPI", false); break; } - case SignalGenerationParameters::SpinEcho: + case SignalGenerationParameters::ConventionalSpinEcho: { PrintToLog("Acquisition type: classic spin echo with cartesian k-space trajectory", false); break; } default: { PrintToLog("Acquisition type: single shot EPI", false); break; } } if(m_Parameters.m_SignalGen.m_tInv>0) PrintToLog("Using inversion pulse with TI " + boost::lexical_cast(m_Parameters.m_SignalGen.m_tInv) + "ms", false); if (m_Parameters.m_SignalGen.m_DoSimulateRelaxation) PrintToLog("Simulating signal relaxation", false); if (m_Parameters.m_SignalGen.m_NoiseVariance>0 && m_Parameters.m_Misc.m_DoAddNoise) PrintToLog("Simulating complex Gaussian noise: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_NoiseVariance), false); if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull() && m_Parameters.m_Misc.m_DoAddDistortions) PrintToLog("Simulating distortions", false); if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging) { if (m_Parameters.m_SignalGen.m_ZeroRinging > 0) PrintToLog("Simulating ringing artifacts by zeroing " + boost::lexical_cast(m_Parameters.m_SignalGen.m_ZeroRinging) + "% of k-space frequencies", false); else PrintToLog("Simulating ringing artifacts by cropping high resolution inputs during k-space simulation", false); } if (m_Parameters.m_Misc.m_DoAddEddyCurrents && m_Parameters.m_SignalGen.m_EddyStrength>0) PrintToLog("Simulating eddy currents: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_EddyStrength), false); if (m_Parameters.m_Misc.m_DoAddSpikes && m_Parameters.m_SignalGen.m_Spikes>0) PrintToLog("Simulating spikes: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_Spikes), false); if (m_Parameters.m_Misc.m_DoAddAliasing && m_Parameters.m_SignalGen.m_CroppingFactor<1.0) PrintToLog("Simulating aliasing: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_CroppingFactor), false); if (m_Parameters.m_Misc.m_DoAddGhosts && m_Parameters.m_SignalGen.m_KspaceLineOffset>0) PrintToLog("Simulating ghosts: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_KspaceLineOffset), false); doubleOutImage = SimulateKspaceAcquisition(m_CompartmentImages); signalScale = 1; // already scaled in SimulateKspaceAcquisition() } else // don't do k-space stuff, just sum compartments { PrintToLog("Summing compartments"); doubleOutImage = m_CompartmentImages.at(0); for (unsigned int i=1; i::New(); adder->SetInput1(doubleOutImage); adder->SetInput2(m_CompartmentImages.at(i)); adder->Update(); doubleOutImage = adder->GetOutput(); } } if (this->GetAbortGenerateData()) { PrintToLog("\n", false, false); PrintToLog("Simulation aborted"); return; } PrintToLog("Finalizing image"); if (m_Parameters.m_SignalGen.m_DoAddDrift && m_Parameters.m_SignalGen.m_Drift>0.0) PrintToLog("Adding signal drift: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_Drift), false); if (signalScale>1) PrintToLog("Scaling signal", false); if (m_Parameters.m_NoiseModel) PrintToLog("Adding noise: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_NoiseVariance), false); 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()); PrintToLog("0% 10 20 30 40 50 60 70 80 90 100%", false, true, false); PrintToLog("|----|----|----|----|----|----|----|----|----|----|\n*", false, false, false); int lastTick = 0; while(!it4.IsAtEnd()) { if (this->GetAbortGenerateData()) { PrintToLog("\n", false, false); PrintToLog("Simulation aborted"); return; } ++disp2; unsigned long newTick = 50*disp2.count()/disp2.expected_count(); for (unsigned long tick = 0; tick<(newTick-lastTick); tick++) PrintToLog("*", false, false, false); lastTick = newTick; typename OutputImageType::IndexType index = it4.GetIndex(); signal = doubleOutImage->GetPixel(index)*signalScale; for (unsigned int i=0; iAddNoise(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); PrintToLog("\n", false); PrintToLog("Finished simulation"); m_TimeProbe.Stop(); if (m_Parameters.m_SignalGen.m_DoAddMotion) { PrintToLog("\nHead motion log:", false); PrintToLog(m_MotionLog, false, false); } if (m_Parameters.m_Misc.m_DoAddSpikes && m_Parameters.m_SignalGen.m_Spikes>0) { PrintToLog("\nSpike log:", false); PrintToLog(m_SpikeLog, false, false); } if (m_Logfile.is_open()) m_Logfile.close(); } template< class PixelType > void TractsToDWIImageFilter< PixelType >::PrintToLog(std::string m, bool addTime, bool linebreak, bool stdOut) { // timestamp if (addTime) { m_Logfile << this->GetTime() << " > "; m_StatusText += this->GetTime() + " > "; if (stdOut) std::cout << this->GetTime() << " > "; } // message if (m_Logfile.is_open()) m_Logfile << m; m_StatusText += m; if (stdOut) std::cout << m; // new line if (linebreak) { if (m_Logfile.is_open()) m_Logfile << "\n"; m_StatusText += "\n"; if (stdOut) std::cout << "\n"; } m_Logfile.flush(); } template< class PixelType > void TractsToDWIImageFilter< PixelType >::SimulateMotion(int g) { if ( m_Parameters.m_SignalGen.m_DoAddMotion && m_Parameters.m_SignalGen.m_DoRandomizeMotion && g>0 && m_Parameters.m_SignalGen.m_MotionVolumes[g-1]) { // The last volume was randomly moved, so we have to reset to fiberbundle and the mask. // Without motion or with linear motion, we keep the last position --> no reset. m_FiberBundleTransformed = m_FiberBundle->GetDeepCopy(); if (m_MaskImageSet) { auto duplicator = itk::ImageDuplicator::New(); duplicator->SetInputImage(m_Parameters.m_SignalGen.m_MaskImage); duplicator->Update(); m_TransformedMaskImage = duplicator->GetOutput(); } } VectorType rotation; VectorType translation; // is motion artifact enabled? // is the current volume g affected by motion? if ( m_Parameters.m_SignalGen.m_DoAddMotion && m_Parameters.m_SignalGen.m_MotionVolumes[g] && g(m_Parameters.m_SignalGen.GetNumVolumes()) ) { // adjust motion transforms if ( m_Parameters.m_SignalGen.m_DoRandomizeMotion ) { // randomly rotation[0] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[0]*2) -m_Parameters.m_SignalGen.m_Rotation[0]; rotation[1] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[1]*2) -m_Parameters.m_SignalGen.m_Rotation[1]; rotation[2] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[2]*2) -m_Parameters.m_SignalGen.m_Rotation[2]; translation[0] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[0]*2) -m_Parameters.m_SignalGen.m_Translation[0]; translation[1] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[1]*2) -m_Parameters.m_SignalGen.m_Translation[1]; translation[2] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[2]*2) -m_Parameters.m_SignalGen.m_Translation[2]; m_FiberBundleTransformed->TransformFibers(rotation[0], rotation[1], rotation[2], translation[0], translation[1], translation[2]); } else { // linearly rotation = m_Parameters.m_SignalGen.m_Rotation / m_NumMotionVolumes; translation = m_Parameters.m_SignalGen.m_Translation / m_NumMotionVolumes; m_MotionCounter++; m_FiberBundleTransformed->TransformFibers(rotation[0], rotation[1], rotation[2], translation[0], translation[1], translation[2]); rotation *= m_MotionCounter; translation *= m_MotionCounter; } MatrixType rotationMatrix = mitk::imv::GetRotationMatrixItk(rotation[0], rotation[1], rotation[2]); MatrixType rotationMatrixInv = mitk::imv::GetRotationMatrixItk(-rotation[0], -rotation[1], -rotation[2]); m_Rotations.push_back(rotationMatrix); m_RotationsInv.push_back(rotationMatrixInv); m_Translations.push_back(translation); // move mask image accoring to new transform 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(); m_TransformedMaskImage->TransformPhysicalPointToIndex(GetMovedPoint(index, true), index); if (m_TransformedMaskImage->GetLargestPossibleRegion().IsInside(index)) m_TransformedMaskImage->SetPixel(index, 100); ++maskIt; } } } else { if (m_Parameters.m_SignalGen.m_DoAddMotion && !m_Parameters.m_SignalGen.m_DoRandomizeMotion && g>0) { rotation = m_Parameters.m_SignalGen.m_Rotation / m_NumMotionVolumes; rotation *= m_MotionCounter; m_Rotations.push_back(m_Rotations.back()); m_RotationsInv.push_back(m_RotationsInv.back()); m_Translations.push_back(m_Translations.back()); } else { rotation.Fill(0.0); VectorType translation; translation.Fill(0.0); MatrixType rotation_matrix; rotation_matrix.SetIdentity(); m_Rotations.push_back(rotation_matrix); m_RotationsInv.push_back(rotation_matrix); m_Translations.push_back(translation); } } if (m_Parameters.m_SignalGen.m_DoAddMotion) { m_MotionLog += boost::lexical_cast(g) + " rotation: " + boost::lexical_cast(rotation[0]) + "," + boost::lexical_cast(rotation[1]) + "," + boost::lexical_cast(rotation[2]) + ";"; m_MotionLog += " translation: " + boost::lexical_cast(m_Translations.back()[0]) + "," + boost::lexical_cast(m_Translations.back()[1]) + "," + boost::lexical_cast(m_Translations.back()[2]) + "\n"; } } template< class PixelType > itk::Point TractsToDWIImageFilter< PixelType >::GetMovedPoint(itk::Index<3>& index, bool forward) { itk::Point transformed_point; float tx = m_Translations.back()[0]; float ty = m_Translations.back()[1]; float tz = m_Translations.back()[2]; if (forward) { m_UpsampledMaskImage->TransformIndexToPhysicalPoint(index, transformed_point); m_FiberBundle->TransformPoint<>(transformed_point, m_Rotations.back(), tx, ty, tz); } else { tx *= -1; ty *= -1; tz *= -1; m_TransformedMaskImage->TransformIndexToPhysicalPoint(index, transformed_point); m_FiberBundle->TransformPoint<>(transformed_point, m_RotationsInv.back(), tx, ty, tz); } return transformed_point; } template< class PixelType > void TractsToDWIImageFilter< PixelType >:: SimulateExtraAxonalSignal(ItkUcharImgType::IndexType& index, itk::Point& volume_fraction_point, double intraAxonalVolume, int g) { int numFiberCompartments = m_Parameters.m_FiberModelList.size(); int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size(); if (m_Parameters.m_SignalGen.m_DoDisablePartialVolume) { // simulate signal for largest non-fiber compartment int max_compartment_index = 0; double max_fraction = 0; if (numNonFiberCompartments>1) { for (int i=0; iSetInputImage(m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()); double compartment_fraction = mitk::imv::GetImageValue(volume_fraction_point, true, m_DoubleInterpolator); if (compartment_fraction<0) mitkThrow() << "Volume fraction image (index " << i << ") contains values less than zero!"; if (compartment_fraction>max_fraction) { max_fraction = compartment_fraction; max_compartment_index = i; } } } DoubleDwiType::Pointer doubleDwi = m_CompartmentImages.at(max_compartment_index+numFiberCompartments); DoubleDwiType::PixelType pix = doubleDwi->GetPixel(index); pix[g] += m_Parameters.m_NonFiberModelList[max_compartment_index]->SimulateMeasurement(g, m_NullDir)*m_VoxelVolume; doubleDwi->SetPixel(index, pix); if (g==0) m_VolumeFractions.at(max_compartment_index+numFiberCompartments)->SetPixel(index, 1); } else { std::vector< double > fractions; if (g==0) m_VolumeFractions.at(0)->SetPixel(index, intraAxonalVolume/m_VoxelVolume); double extraAxonalVolume = m_VoxelVolume-intraAxonalVolume; // non-fiber volume if (extraAxonalVolume<0) { if (extraAxonalVolume<-0.001) MITK_ERROR << "Corrupted intra-axonal signal voxel detected. Fiber volume larger voxel volume! " << m_VoxelVolume << "<" << intraAxonalVolume; extraAxonalVolume = 0; } double interAxonalVolume = 0; if (numFiberCompartments>1) interAxonalVolume = extraAxonalVolume * intraAxonalVolume/m_VoxelVolume; // inter-axonal fraction of non fiber compartment double nonFiberVolume = extraAxonalVolume - interAxonalVolume; // rest of compartment if (nonFiberVolume<0) { if (nonFiberVolume<-0.001) MITK_ERROR << "Corrupted signal voxel detected. Fiber volume larger voxel volume!"; nonFiberVolume = 0; interAxonalVolume = extraAxonalVolume; } double compartmentSum = intraAxonalVolume; fractions.push_back(intraAxonalVolume/m_VoxelVolume); // rescale extra-axonal fiber signal for (int i=1; iGetVolumeFractionImage()!=nullptr) { m_DoubleInterpolator->SetInputImage(m_Parameters.m_FiberModelList[i]->GetVolumeFractionImage()); interAxonalVolume = mitk::imv::GetImageValue(volume_fraction_point, true, m_DoubleInterpolator)*m_VoxelVolume; if (interAxonalVolume<0) mitkThrow() << "Volume fraction image (index " << i+1 << ") contains negative values!"; } DoubleDwiType::PixelType pix = m_CompartmentImages.at(i)->GetPixel(index); pix[g] *= interAxonalVolume; m_CompartmentImages.at(i)->SetPixel(index, pix); compartmentSum += interAxonalVolume; fractions.push_back(interAxonalVolume/m_VoxelVolume); if (g==0) m_VolumeFractions.at(i)->SetPixel(index, interAxonalVolume/m_VoxelVolume); } for (int i=0; iGetVolumeFractionImage()!=nullptr) { m_DoubleInterpolator->SetInputImage(m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()); volume = mitk::imv::GetImageValue(volume_fraction_point, true, m_DoubleInterpolator)*m_VoxelVolume; if (volume<0) mitkThrow() << "Volume fraction image (index " << numFiberCompartments+i+1 << ") contains negative values (non-fiber compartment)!"; if (m_UseRelativeNonFiberVolumeFractions) volume *= nonFiberVolume/m_VoxelVolume; } DoubleDwiType::PixelType pix = m_CompartmentImages.at(i+numFiberCompartments)->GetPixel(index); pix[g] += m_Parameters.m_NonFiberModelList[i]->SimulateMeasurement(g, m_NullDir)*volume; m_CompartmentImages.at(i+numFiberCompartments)->SetPixel(index, pix); compartmentSum += volume; fractions.push_back(volume/m_VoxelVolume); if (g==0) m_VolumeFractions.at(i+numFiberCompartments)->SetPixel(index, volume/m_VoxelVolume); } if (compartmentSum/m_VoxelVolume>1.05) { MITK_ERROR << "Compartments do not sum to 1 in voxel " << index << " (" << compartmentSum/m_VoxelVolume << ")"; for (auto val : fractions) MITK_ERROR << val; } } } template< class PixelType > itk::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/IODataStructures/mitkFiberfoxParameters.h b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h index bcb285d2f4..ec61f3ea82 100644 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h @@ -1,333 +1,333 @@ /*=================================================================== 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 #include #include #include #include #include namespace mitk { class MITKFIBERTRACKING_EXPORT FiberfoxParameters; /** Signal generation */ class MITKFIBERTRACKING_EXPORT SignalGenerationParameters { friend FiberfoxParameters; public: typedef itk::Image ItkFloatImgType; typedef itk::Image ItkUcharImgType; typedef itk::Vector GradientType; typedef std::vector GradientListType; enum CoilSensitivityProfile : int { COIL_CONSTANT, COIL_LINEAR, COIL_EXPONENTIAL }; enum AcquisitionType : int { SingleShotEpi, - SpinEcho + ConventionalSpinEcho }; SignalGenerationParameters() : m_AcquisitionType(SignalGenerationParameters::SingleShotEpi) , m_SignalScale(100) , m_tEcho(100) , m_tRep(4000) , m_tInv(0) , m_tLine(1) , m_tInhom(50) , m_ReversePhase(false) , m_PartialFourier(1.0) , m_NoiseVariance(0.001) , m_NumberOfCoils(1) , m_CoilSensitivityProfile(SignalGenerationParameters::COIL_CONSTANT) , m_SimulateKspaceAcquisition(false) , m_AxonRadius(0) , m_DoDisablePartialVolume(false) , m_Spikes(0) , m_SpikeAmplitude(1) , m_KspaceLineOffset(0) , m_EddyStrength(300) , m_Tau(70) , m_CroppingFactor(1) , m_Drift(0.06) , m_DoAddGibbsRinging(false) , m_ZeroRinging(0) , m_DoSimulateRelaxation(true) , m_DoAddMotion(false) , m_DoRandomizeMotion(true) , m_DoAddDrift(false) , m_FrequencyMap(nullptr) , m_MaskImage(nullptr) , m_Bvalue(1000) { 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_CroppedRegion; ///< Image size with reduced FOV. 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 */ AcquisitionType m_AcquisitionType; ///< determines k-space trajectory and maximum echo position(s) float m_SignalScale; ///< Scaling factor for output signal (before noise is added). float m_tEcho; ///< Echo time TE. float m_tRep; ///< Echo time TR. float m_tInv; ///< Inversion time float m_tLine; ///< k-space line readout time (dwell time). float m_tInhom; ///< T2' bool m_ReversePhase; ///< If true, the phase readout direction will be inverted (-y instead of y) float m_PartialFourier; ///< Partial fourier factor (0.5-1) float m_NoiseVariance; ///< Variance of complex gaussian noise unsigned int m_NumberOfCoils; ///< Number of coils in multi-coil acquisition CoilSensitivityProfile m_CoilSensitivityProfile; ///< Choose between constant, linear or exponential sensitivity profile of the used coils 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. /** Artifacts and other effects */ unsigned int m_Spikes; ///< Number of spikes randomly appearing in the image float m_SpikeAmplitude; ///< amplitude of spikes relative to the largest signal intensity (magnitude of complex) float m_KspaceLineOffset; ///< Causes N/2 ghosts. Larger offset means stronger ghost. float m_EddyStrength; ///< Strength of eddy current induced gradients in mT/m. float m_Tau; ///< Eddy current decay constant (in ms) float m_CroppingFactor; ///< FOV size in y-direction is multiplied by this factor. Causes aliasing artifacts. float m_Drift; ///< Global signal decrease by the end of the acquisition. bool m_DoAddGibbsRinging; ///< Add Gibbs ringing artifact int m_ZeroRinging; ///< If > 0, ringing is simulated by by setting the defined percentage of higher frequencies to 0 in k-space. Otherwise, the input to the k-space simulation is generated with twice the resolution and cropped during k-space simulation (much slower). bool m_DoSimulateRelaxation; ///< Add T2 relaxation effects bool m_DoAddMotion; ///< Enable motion artifacts. bool m_DoRandomizeMotion; ///< Toggles between random and linear motion. bool m_DoAddDrift; ///< Add quadratic signal drift. std::vector< bool > m_MotionVolumes; ///< Indicates the image volumes that are affected by motion ///< with positive numbers, inverted logic with negative numbers. itk::Vector m_Translation; ///< Maximum translational motion. itk::Vector m_Rotation; ///< Maximum rotational motion. ItkFloatImgType::Pointer m_FrequencyMap; ///< If != nullptr, distortions are added to the image using this frequency map. ItkUcharImgType::Pointer m_MaskImage; ///< Signal is only genrated inside of the mask image. std::vector< int > GetBaselineIndices(); ///< Returns list of nun-diffusion-weighted image volume indices unsigned int GetFirstBaselineIndex(); ///< Returns index of first non-diffusion-weighted image volume bool IsBaselineIndex(unsigned int idx); ///< Checks if image volume with given index is non-diffusion-weighted volume or not. unsigned int GetNumWeightedVolumes(); ///< Get number of diffusion-weighted image volumes unsigned int GetNumBaselineVolumes(); ///< Get number of non-diffusion-weighted image volumes unsigned int GetNumVolumes(); ///< Get number of baseline and diffusion-weighted image volumes GradientListType GetGradientDirections(); ///< Return gradient direction container mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer GetItkGradientContainer(); GradientType GetGradientDirection(unsigned int i); std::vector< int > GetBvalues(); ///< Returns a vector with all unique b-values (determined by the gradient magnitudes) double GetBvalue(); void ApplyDirectionMatrix(); 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. double m_Bvalue; ///< Acquisition b-value void SetNumWeightedVolumes(int numGradients); ///< Automaticall calls GenerateGradientHalfShell() afterwards. void SetGradienDirections(GradientListType gradientList); void SetGradienDirections(mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradientList); void GenerateGradientHalfShell(); ///< Generates half shell of gradient directions (with m_NumGradients non-zero directions) }; /** Fiber generation */ class MITKFIBERTRACKING_EXPORT FiberGenerationParameters { public: enum FiberDistribution { DISTRIBUTE_UNIFORM, // distribute fibers uniformly in the ROIs DISTRIBUTE_GAUSSIAN // distribute fibers using a 2D gaussian }; typedef std::vector< std::vector< mitk::PlanarEllipse::Pointer > > FiducialListType; typedef std::vector< std::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 MITKFIBERTRACKING_EXPORT MiscFiberfoxParameters { public: MiscFiberfoxParameters() : m_ResultNode(DataNode::New()) , m_ParentNode(nullptr) , m_SignalModelString("") , m_ArtifactModelString("") , m_OutputPath("/tmp/") , m_OutputPrefix("fiberfox") , m_AfterSimulationMessage("") , m_BvalsFile("") , m_BvecsFile("") , m_CheckOutputVolumeFractionsBox(false) , m_CheckAdvancedSignalOptionsBox(false) , m_DoAddNoise(false) , m_DoAddGhosts(false) , m_DoAddAliasing(false) , m_DoAddSpikes(false) , m_DoAddEddyCurrents(false) , m_DoAddDistortions(false) , m_MotionVolumesBox("random") , 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. std::string m_SignalModelString; ///< Appendet to the name of the result node std::string m_ArtifactModelString; ///< Appendet to the name of the result node std::string m_OutputPath; ///< Image is automatically saved to the specified folder after simulation is finished. std::string m_OutputPrefix; /** Prefix for filename of output files and logfile. */ std::string m_AfterSimulationMessage; ///< Store messages that are displayed after the simulation has finished (e.g. warnings, automatic parameter adjustments etc.) std::string m_BvalsFile; std::string m_BvecsFile; /** member variables that store the check-state of GUI checkboxes */ // image generation bool m_CheckOutputVolumeFractionsBox; bool m_CheckAdvancedSignalOptionsBox; bool m_DoAddNoise; bool m_DoAddGhosts; bool m_DoAddAliasing; bool m_DoAddSpikes; bool m_DoAddEddyCurrents; bool m_DoAddDistortions; std::string m_MotionVolumesBox; // fiber generation bool m_CheckRealTimeFibersBox; bool m_CheckAdvancedFiberOptionsBox; bool m_CheckConstantRadiusBox; bool m_CheckIncludeFiducialsBox; }; /** * \brief Datastructure to manage the Fiberfox signal generation parameters. * */ class MITKFIBERTRACKING_EXPORT FiberfoxParameters { public: typedef itk::Image ItkFloatImgType; typedef itk::Image ItkDoubleImgType; typedef itk::Image ItkUcharImgType; typedef DiffusionSignalModel DiffusionModelType; typedef std::vector< DiffusionModelType* > DiffusionModelListType; typedef DiffusionNoiseModel NoiseModelType; FiberfoxParameters(); FiberfoxParameters(const FiberfoxParameters ¶ms); ~FiberfoxParameters(); /** Not templated parameters */ FiberGenerationParameters m_FiberGen; ///< Fiber generation parameters SignalGenerationParameters m_SignalGen; ///< Signal generation parameters MiscFiberfoxParameters m_Misc; ///< GUI realted and I/O parameters /** Templated parameters */ DiffusionModelListType m_FiberModelList; ///< Intra- and inter-axonal compartments. DiffusionModelListType m_NonFiberModelList; ///< Extra-axonal compartments. std::shared_ptr< NoiseModelType > m_NoiseModel; ///< If != nullptr, noise is added to the image. void GenerateGradientHalfShell(); void SetNumWeightedVolumes(int numGradients); ///< Automaticall calls GenerateGradientHalfShell() afterwards. void SetGradienDirections(mitk::SignalGenerationParameters::GradientListType gradientList); void SetGradienDirections(mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradientList); void SetBvalue(double Bvalue); void UpdateSignalModels(); void ClearFiberParameters(); void ClearSignalParameters(); void ApplyDirectionMatrix(); void PrintSelf(); ///< Print parameters to stdout. void SaveParameters(std::string filename); ///< Save image generation parameters to .ffp file. void LoadParameters(std::string filename, bool fix_seed=false); ///< Load image generation parameters from .ffp file. template< class ParameterType > ParameterType ReadVal(boost::property_tree::ptree::value_type const& v, std::string tag, ParameterType defaultValue, bool essential=false); std::string m_MissingTags; }; } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/files.cmake b/Modules/DiffusionImaging/FiberTracking/files.cmake index a93ec3dcd5..5567d247a5 100644 --- a/Modules/DiffusionImaging/FiberTracking/files.cmake +++ b/Modules/DiffusionImaging/FiberTracking/files.cmake @@ -1,108 +1,108 @@ set(CPP_FILES mitkFiberTrackingModuleActivator.cpp ## IO datastructures IODataStructures/FiberBundle/mitkFiberBundle.cpp IODataStructures/FiberBundle/mitkTrackvis.cpp IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp IODataStructures/mitkTractographyForest.cpp IODataStructures/mitkFiberfoxParameters.cpp # Interactions # Tractography Algorithms/GibbsTracking/mitkParticleGrid.cpp Algorithms/GibbsTracking/mitkMetropolisHastingsSampler.cpp Algorithms/GibbsTracking/mitkEnergyComputer.cpp Algorithms/GibbsTracking/mitkGibbsEnergyComputer.cpp Algorithms/GibbsTracking/mitkFiberBuilder.cpp Algorithms/GibbsTracking/mitkSphereInterpolator.cpp Algorithms/itkStreamlineTrackingFilter.cpp Algorithms/TrackingHandlers/mitkTrackingDataHandler.cpp Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.cpp ) set(H_FILES # DataStructures -> FiberBundle IODataStructures/FiberBundle/mitkFiberBundle.h IODataStructures/FiberBundle/mitkTrackvis.h IODataStructures/mitkFiberfoxParameters.h IODataStructures/mitkTractographyForest.h # Algorithms Algorithms/itkTractDensityImageFilter.h Algorithms/itkTractsToFiberEndingsImageFilter.h Algorithms/itkTractsToRgbaImageFilter.h Algorithms/itkTractsToVectorImageFilter.h Algorithms/itkEvaluateDirectionImagesFilter.h Algorithms/itkEvaluateTractogramDirectionsFilter.h Algorithms/itkFiberCurvatureFilter.h Algorithms/itkFitFibersToImageFilter.h Algorithms/itkTractClusteringFilter.h Algorithms/itkTractDistanceFilter.h Algorithms/itkFiberExtractionFilter.h Algorithms/itkTdiToVolumeFractionFilter.h # Tractography Algorithms/TrackingHandlers/mitkTrackingDataHandler.h Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.h Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.h Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.h Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.h Algorithms/itkGibbsTrackingFilter.h Algorithms/itkStochasticTractographyFilter.h Algorithms/GibbsTracking/mitkParticle.h Algorithms/GibbsTracking/mitkParticleGrid.h Algorithms/GibbsTracking/mitkMetropolisHastingsSampler.h Algorithms/GibbsTracking/mitkSimpSamp.h Algorithms/GibbsTracking/mitkEnergyComputer.h Algorithms/GibbsTracking/mitkGibbsEnergyComputer.h Algorithms/GibbsTracking/mitkSphereInterpolator.h Algorithms/GibbsTracking/mitkFiberBuilder.h Algorithms/itkStreamlineTrackingFilter.h # Clustering Algorithms/ClusteringMetrics/mitkClusteringMetric.h Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanMean.h Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanMax.h Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanStd.h Algorithms/ClusteringMetrics/mitkClusteringMetricAnatomic.h Algorithms/ClusteringMetrics/mitkClusteringMetricScalarMap.h Algorithms/ClusteringMetrics/mitkClusteringMetricInnerAngles.h Algorithms/ClusteringMetrics/mitkClusteringMetricLength.h # Fiberfox Fiberfox/itkFibersFromPlanarFiguresFilter.h Fiberfox/itkTractsToDWIImageFilter.h Fiberfox/itkKspaceImageFilter.h Fiberfox/itkDftImageFilter.h Fiberfox/itkFieldmapGeneratorFilter.h Fiberfox/itkRandomPhantomFilter.h Fiberfox/SignalModels/mitkDiffusionSignalModel.h Fiberfox/SignalModels/mitkTensorModel.h Fiberfox/SignalModels/mitkBallModel.h Fiberfox/SignalModels/mitkDotModel.h Fiberfox/SignalModels/mitkAstroStickModel.h Fiberfox/SignalModels/mitkStickModel.h Fiberfox/SignalModels/mitkRawShModel.h Fiberfox/SignalModels/mitkDiffusionNoiseModel.h Fiberfox/SignalModels/mitkRicianNoiseModel.h Fiberfox/SignalModels/mitkChiSquareNoiseModel.h Fiberfox/Sequences/mitkAcquisitionType.h Fiberfox/Sequences/mitkSingleShotEpi.h - Fiberfox/Sequences/mitkCartesianReadout.h + Fiberfox/Sequences/mitkConventionalSpinEcho.h ) set(RESOURCE_FILES # Binary directory resources FiberTrackingLUTBaryCoords.bin FiberTrackingLUTIndices.bin ) diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberfox/src/internal/QmitkFiberfoxViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberfox/src/internal/QmitkFiberfoxViewControls.ui index 3966ac3473..a3373729d4 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberfox/src/internal/QmitkFiberfoxViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberfox/src/internal/QmitkFiberfoxViewControls.ui @@ -1,2793 +1,2845 @@ QmitkFiberfoxViewControls 0 0 - 463 - 2744 + 490 + 2775 Form - - + + - - - - Save Parameters - - - - :/QmitkDiffusionImaging/general_icons/download.ico:/QmitkDiffusionImaging/general_icons/download.ico - - - - - - - QFrame::NoFrame + QGroupBox { + background-color: transparent; +} - - QFrame::Raised + + Intra-axonal Compartment - + - 0 + 6 - 0 + 6 - 0 + 6 - 0 + 6 + + + - - - 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 + Select signal model for intra-axonal compartment. + + + Stick Model + + + + + Zeppelin Model + + + + + Tensor Model + + + + + Prototype Signal + + - - - - QGroupBox { - background-color: transparent; -} + + + + + + + QFrame::NoFrame - - Input Data + + QFrame::Raised - + - 6 + 0 - 6 + 0 - 6 + 0 - 6 + 0 - - - - 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 + Volume Fraction: - - + + - <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 + 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 - - - - - - - - Courier - 7 - - - - true - - - - - + + QGroupBox { background-color: transparent; } - Intra-axonal Compartment + Inter-axonal Compartment - + 6 6 6 6 - - - - - - - - - - + Select signal model for intra-axonal compartment. - Stick Model + -- - Zeppelin Model + Stick Model - Tensor Model + Zeppelin Model - Prototype Signal + 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. - - - - - - - Load Parameters - - - - :/QmitkDiffusionImaging/general_icons/upload.ico:/QmitkDiffusionImaging/general_icons/upload.ico - - - - - + + QGroupBox { background-color: transparent; } - Extra-axonal Compartments + Image Settings - + 6 6 6 6 - - + + + + color: rgb(255, 0, 0); + + + Using geometry of selected image! + + + + + QFrame::NoFrame QFrame::Raised - + 0 0 0 0 - - - - Volume Fraction: + + 6 + + + + + TE in milliseconds + + + 1 + + + 10000 + + + 1 + + + 100 - - - - + + + + Single Shot EPI + + + + + Conventional Spin Echo + + - - - - - - - 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: + Fiber Radius: - - + + - Optional! If no volume fraction map for this compartment is set, the corresponding volume fractions are calculated from the input fibers. + Number of coil elements used for the acquisiton. + + + 1 + + + 128 + + + 1 + + + 1 + + + + + + + TE in milliseconds + + + 1 + + + 999999999 + + + 1 + + + 100 + + + + + + + Output phase image and volume fraction maps. + + + Output Additional Images + + + false + + + + + + + Dwell time (time to read one line in k-space) in ms. + + + 100.000000000000000 + + + 0.100000000000000 + + + 1.000000000000000 + + + + + + + Partial fourier factor (0.5-1) + + + 3 + + + 0.500000000000000 + + + 1.000000000000000 + + + 0.100000000000000 + + + 1.000000000000000 - - - - - - - - - - - - - - - - Select signal model for extra-axonal compartment. - - - - -- - - - - - Ball Model - - - - - Astrosticks Model - - - - - Dot Model - - - - - Prototype Signal - - - - - - - - - - - - - - QGroupBox { - background-color: transparent; -} - - - Noise and other Artifacts - - - - 6 - - - 6 - - - 6 - - - 6 - - - - - Qt::Horizontal - - - - - - - true - - - QFrame::NoFrame - - - QFrame::Raised - - - - 0 - - - 6 - - - 0 - - - 0 - - - 6 - - + + + Acquisition Type: + + + + + - Toggle between random movement and linear movement. + Fiber radius used to calculate volume fractions (in µm). Set to 0 for automatic radius estimation. + + + 9999.000000000000000 + + + + + + + + + + + + + - Randomize motion + Partial Fourier: + + + false + + + + + + + Output one image per compartment containing the corresponding volume fractions per voxel. + + + Reverse Phase Encoding Direction - true + false - - - QGroupBox { - background-color: transparent; -} + + + - - Rotation + + + + + + + + <html><head/><body><p>Coil Sensitivity:</p></body></html> + + + false + + + + + + + + + + + + + + + + <html><head/><body><p>Number of Channels:</p></body></html> + + + false + + + + + + + + + + + + + + + + <html><head/><body><p>Repetition Time <span style=" font-style:italic;">TR</span>: </p></body></html> + + + false + + + + + + + TR in milliseconds + + + 1 + + + 999999999 + + + 1 + + + 4000 + + + + + + + Signal Scale: + + + + + + + Relaxation time due to magnetic field inhomogeneities (T2', in milliseconds). + + + 1 + + + 10000 + + + 1 + + + 50 + + + + + + + + + + + + + + + + <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 + + + + + + + <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 + + + true + + + + + + + + + + + + + + + + <html><head/><body><p>Echo Time <span style=" font-style:italic;">TE</span>: </p></body></html> + + + false + + + + + + + + + + + + + + + + Dwell Time: + + + false + + + + + + + + Constant + + + + + Linear + + + + + Exponential + + + + + + + + Disable partial volume. Treat voxel content as fiber-only if at least one fiber is present. + + + Disable Partial Volume Effects + + + false + + + + + + + + + + + + + + + + <html><head/><body><p>Inversion Time <span style=" font-style:italic;">TI</span>: </p></body></html> + + + false + + + + + + + Inversion time (in ms) for inversion recovery sequences. If 0, no inversion pulse is simulated. + + + 0 + + + 999999999 + + + 1 + + + 0 + + + + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + 0 + + + 0 + + + 0 + + + 6 + + + + + + + + + + + + + + <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 - - - 6 - - - 9 - - - 6 - - - 6 - - - - - - - - - - - - - - Degree: - - - false - - - - - - - - - - - - - - - - x - - - false - - - - - - - - - - - - - - - - Axis: - - - false - - - - - - - Maximum rotation around x-axis. - - - 1 - - - -360.000000000000000 - - - 360.000000000000000 - - - 1.000000000000000 - - - 0.000000000000000 - - - - - - - Maximum rotation around z-axis. - - - 1 - - - -360.000000000000000 - - - 360.000000000000000 - - - 1.000000000000000 - - - 15.000000000000000 - - - - - - - - - - - - - - - - y - - - false - - - - - - - - - - - - - - - - z - - - false - - - - - - - Maximum rotation around y-axis. - - - 1 - - - -360.000000000000000 - - - 360.000000000000000 - - - 1.000000000000000 - - - 0.000000000000000 - - - - - - - - QGroupBox { - background-color: transparent; -} + + + + Gradient Directions: - - Translation + + + + + + Number of gradient directions distributed over the half sphere. + + + 0 + + + 10000 + + + 1 + + + 30 + + + + + + + + + + Advanced Options + + + + + + + color: rgb(255, 0, 0); + + + Using gradients of selected DWI! + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + 0 + + + 0 + + + 0 + + + + + 3 + + + 0.100000000000000 + + + 50.000000000000000 + + + 0.100000000000000 + + + 2.000000000000000 - - - 6 - - - 6 - - - 6 - - - - - - - - - - - - - - Distance: - - - false - - - - - - - - - - - - - - - - x - - - false - - - - - - - - - - - - - - - - y - - - false - - - - - - - - - - - - - - - - Axis: - - - false - - - - - - - - - - - - - - - - z - - - false - - - - - - - Maximum translation along x-axis. - - - 1 - - - -1000.000000000000000 - - - 1000.000000000000000 - - - 1.000000000000000 - - - 0.000000000000000 - - - - - - - Maximum translation along y-axis. - - - 1 - - - -1000.000000000000000 - - - 1000.000000000000000 - - - 1.000000000000000 - - - 0.000000000000000 - - - - - - - Maximum translation along z-axis. - - - 1 - - - -1000.000000000000000 - - - 1000.000000000000000 - - - 1.000000000000000 - - - 0.000000000000000 - - - - - - - QFrame::NoFrame + + + Image Spacing: - - QFrame::Raised + + + + + + 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 + + + 20 + + + + + + + Fiber sampling factor which determines the accuracy of the calculated fiber and non-fiber volume fractions. + + + 1 + + + 1000 + + + 1 + + + 20 + + + + + + + Fiber sampling factor which determines the accuracy of the calculated fiber and non-fiber volume fractions. + + + 1 + + + 1000 + + + 1 + + + 3 - - - 0 - - - 0 - - - 0 - - - 0 - - - - - Motion volumes: - - - - - - - Type in the volume indices that should be affected by motion (e.g. "0 3 7" whithout quotation marks). Leave blank for motion in all volumes. Type in "random" to randomly select volumes for motion. A list of negative numbers (e.g. -1 -2 -3) excludes volumes (e.g. 1 2 3) selects all remaining volumes. - - - random - - - - - - + + + + Use bvals/bvecs files + + + + + QFrame::NoFrame QFrame::Raised - + 0 0 0 0 - - + + - Num. Spikes: + ... + + + + + + + ... - + + + - + + + false + + + + + - The number of randomly occurring signal spikes. + - - 1 + + + + + + + + Bvecs: + + + false - - + + - Spike amplitude relative to the largest signal amplitude of the corresponding k-space slice. + - - 0.100000000000000 + + - - 0.100000000000000 + + + + + Bvals: + + + false - - + + - Scale: + - + + + false - - - - true - + + + + + + + QGroupBox { + background-color: transparent; +} + + + Extra-axonal Compartments + + + + 6 + + + 6 + + + 6 + + + 6 + + + QFrame::NoFrame QFrame::Raised - - - 6 - + 0 0 0 0 - - - - - - - - - - + - Shrink FOV (%): - - - false + Volume Fraction: - + - Shrink FOV by this percentage. - - - 1 - - - 0.000000000000000 - - - 90.000000000000000 - - - 0.100000000000000 - - - 25.000000000000000 + - - - - Qt::Horizontal + + + + Select signal model for extra-axonal compartment. + + + Ball Model + + + + + Astrosticks Model + + + + + Dot Model + + + + + Prototype Signal + + + + + + + + - - - true + + + + + + + + + Qt::Horizontal + + + + QFrame::NoFrame QFrame::Raised - - - 6 - + 0 0 0 0 - - - - - - - - - - + - Signal Reduction (%): - - - false + Volume Fraction: - + - Global signal in last simulated volume is specified percentage lower than in the first volume. - - - 1 - - - 100.000000000000000 - - - 1.000000000000000 - - - 6.000000000000000 + Optional! If no volume fraction map for this compartment is set, the corresponding volume fractions are calculated from the input fibers. + + + + + + + + + - - - true - - - QFrame::NoFrame - - - QFrame::Raised + + + Select signal model for extra-axonal compartment. - - - 6 + + + -- - - 0 + + + + Ball Model - - 0 + + + + Astrosticks Model - - 0 + + + + Dot Model - - 0 + + + + Prototype Signal - - - - - - - - - - - - - Frequency Map: - - - false - - - - - - - Select image specifying the frequency inhomogeneities (in Hz). - - - - + - - + + + + + + + + + + QGroupBox { + background-color: transparent; +} + + + Noise and other Artifacts + + + + 6 + + + 6 + + + 6 + + + 6 + + + + + Qt::Horizontal + + + + + true QFrame::NoFrame QFrame::Raised - - - QFormLayout::AllNonFixedFieldsGrow - - - 6 - + 0 - 0 + 6 0 0 + + 6 + - + - - - - - - - + Toggle between random movement and linear movement. - Gradient: + Randomize motion - - false + + true - - - - Eddy current induced magnetic field gradient (in mT/m). - - - 4 - - - 1000.000000000000000 - - - 0.001000000000000 - - - 0.010000000000000 + + + + QGroupBox { + background-color: transparent; +} - - - - - - - - - Qt::Horizontal - - - - - - - - - - Add Eddy Current Effects - - - false - - - - - - - Add Distortions - - - false - - - - - - - Add Spikes - - - false - - - - - - - Add Signal Drift - - - false - - - - - - - QFrame::NoFrame - - - QFrame::Raised - - - - 0 - - - 0 - - - 0 - - - 0 - - - - - Variance: + + Rotation + + + 6 + + + 9 + + + 6 + + + 6 + + + + + + + + + + + + + + Degree: + + + false + + + + + + + + + + + + + + + + x + + + false + + + + + + + + + + + + + + + + Axis: + + + false + + + + + + + Maximum rotation around x-axis. + + + 1 + + + -360.000000000000000 + + + 360.000000000000000 + + + 1.000000000000000 + + + 0.000000000000000 + + + + + + + Maximum rotation around z-axis. + + + 1 + + + -360.000000000000000 + + + 360.000000000000000 + + + 1.000000000000000 + + + 15.000000000000000 + + + + + + + + + + + + + + + + y + + + false + + + + + + + + + + + + + + + + z + + + false + + + + + + + Maximum rotation around y-axis. + + + 1 + + + -360.000000000000000 + + + 360.000000000000000 + + + 1.000000000000000 + + + 0.000000000000000 + + + + - - - - Variance of selected noise distribution. - - - 10 - - - 0.000000000000000 - - - 999999999.000000000000000 - - - 0.001000000000000 - - - 50.000000000000000 + + + + QGroupBox { + background-color: transparent; +} - - - - - - Distribution: + + Translation + + + 6 + + + 6 + + + 6 + + + + + + + + + + + + + + Distance: + + + false + + + + + + + + + + + + + + + + x + + + false + + + + + + + + + + + + + + + + y + + + false + + + + + + + + + + + + + + + + Axis: + + + false + + + + + + + + + + + + + + + + z + + + false + + + + + + + Maximum translation along x-axis. + + + 1 + + + -1000.000000000000000 + + + 1000.000000000000000 + + + 1.000000000000000 + + + 0.000000000000000 + + + + + + + Maximum translation along y-axis. + + + 1 + + + -1000.000000000000000 + + + 1000.000000000000000 + + + 1.000000000000000 + + + 0.000000000000000 + + + + + + + Maximum translation along z-axis. + + + 1 + + + -1000.000000000000000 + + + 1000.000000000000000 + + + 1.000000000000000 + + + 0.000000000000000 + + + + - - - - Noise distribution + + + + QFrame::NoFrame - - - Complex Gaussian + + QFrame::Raised + + + + 0 - - - - Rician + + 0 - - - - - - - - - - Qt::Horizontal - - - - - - - Qt::Horizontal - - - - - - - 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 - - - - - - - - - - Add Motion Artifacts - - - false - - - - - - - Add N/2 Ghosts - - - false - - - - - - - Qt::Horizontal - - - - - - - Add ringing artifacts occuring at strong edges in the image. - - - Add Gibbs Ringing - - - false - - - - - - - Qt::Horizontal - - - - - - - Add Noise - - - false - - - - - - - Qt::Horizontal - - - - - - - Add Aliasing - - - false - - - - - - - If > 0, ringing is simulated by by setting the defined percentage of higher frequencies to 0 in k-space. Otherwise, the input to the k-space simulation is generated with twice the resolution and cropped during k-space simulation (much slower). - - - 100 - - - 10 - - - - - - - - - - QGroupBox { - background-color: transparent; -} - - - Image Settings - - - - 6 - - - 6 - - - 6 - - - 6 - - - - - color: rgb(255, 0, 0); - - - Using geometry of selected image! - + + 0 + + + 0 + + + + + Motion volumes: + + + + + + + Type in the volume indices that should be affected by motion (e.g. "0 3 7" whithout quotation marks). Leave blank for motion in all volumes. Type in "random" to randomly select volumes for motion. A list of negative numbers (e.g. -1 -2 -3) excludes volumes (e.g. 1 2 3) selects all remaining volumes. + + + random + + + + + + + - - + + QFrame::NoFrame QFrame::Raised - + 0 0 0 0 - - 6 - - - - - TE in milliseconds - - - 1 - - - 10000 - - - 1 - - - 100 - - - - - - - - Single Shot EPI - - - - - Spin Echo - - - - - - + + - Fiber Radius: + Num. Spikes: - - + + - Number of coil elements used for the acquisiton. - - - 1 - - - 128 - - - 1 + The number of randomly occurring signal spikes. 1 - - - - TE in milliseconds - - - 1 - - - 999999999 - - - 1 - - - 100 - - - - - - - Output phase image and volume fraction maps. - - - Output Additional Images - - - false - - - - - + + - Dwell time (time to read one line in k-space) in ms. - - - 100.000000000000000 + Spike amplitude relative to the largest signal amplitude of the corresponding k-space slice. 0.100000000000000 - 1.000000000000000 - - - - - - - Partial fourier factor (0.5-1) - - - 3 - - - 0.500000000000000 - - - 1.000000000000000 - - 0.100000000000000 - - 1.000000000000000 - - - + + - Acquisition Type: - - - - - - - Fiber radius used to calculate volume fractions (in µm). Set to 0 for automatic radius estimation. - - - 9999.000000000000000 + Scale: - - + + + + + + + true + + + QFrame::NoFrame + + + QFrame::Raised + + + + 6 + + + 0 + + + 0 + + + 0 + + + 0 + + + - Partial Fourier: + Shrink FOV (%): false - - - - Output one image per compartment containing the corresponding volume fractions per voxel. - - - Reverse Phase Encoding Direction - - - false - - - - - + + - - - - - - - - - - <html><head/><body><p>Coil Sensitivity:</p></body></html> - - - false + Shrink FOV by this percentage. - - - - - - + + 1 - - + + 0.000000000000000 - - + + 90.000000000000000 - - <html><head/><body><p>Number of Channels:</p></body></html> + + 0.100000000000000 - - false + + 25.000000000000000 - - + + + + + + + Qt::Horizontal + + + + + + + true + + + QFrame::NoFrame + + + QFrame::Raised + + + + 6 + + + 0 + + + 0 + + + 0 + + + 0 + + + - <html><head/><body><p>Repetition Time <span style=" font-style:italic;">TR</span>: </p></body></html> + Signal Reduction (%): false - - - - TR in milliseconds - - - 1 - - - 999999999 - - - 1 - - - 4000 - - - - - - - Signal Scale: - - - - - + + - Relaxation time due to magnetic field inhomogeneities (T2', in milliseconds). + Global signal in last simulated volume is specified percentage lower than in the first volume. - + 1 - 10000 + 100.000000000000000 - 1 + 1.000000000000000 - 50 + 6.000000000000000 - - + + + + + + + true + + + QFrame::NoFrame + + + QFrame::Raised + + + + 6 + + + 0 + + + 0 + + + 0 + + + 0 + + + - <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> + Frequency Map: false - - + + - <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 - - - true + Select image specifying the frequency inhomogeneities (in Hz). - - + + + + + + + true + + + QFrame::NoFrame + + + QFrame::Raised + + + + QFormLayout::AllNonFixedFieldsGrow + + + 6 + + + 0 + + + 0 + + + 0 + + + 0 + + + - <html><head/><body><p>Echo Time <span style=" font-style:italic;">TE</span>: </p></body></html> + Gradient: false - - + + - + Eddy current induced magnetic field gradient (in mT/m). - - + + 4 - - + + 1000.000000000000000 - - Dwell Time: + + 0.001000000000000 - - false + + 0.010000000000000 - - - - - Constant - - - - - Linear - - - - - Exponential - - - - - - - - Disable partial volume. Treat voxel content as fiber-only if at least one fiber is present. - + + + + + + + Qt::Horizontal + + + + + + + + + + Add Eddy Current Effects + + + false + + + + + + + Add Distortions + + + false + + + + + + + Add Spikes + + + false + + + + + + + Add Signal Drift + + + false + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + 0 + + + 0 + + + 0 + + + - Disable Partial Volume Effects - - - false + Variance: - - + + - - - - - - - - - - <html><head/><body><p>Inversion Time <span style=" font-style:italic;">TI</span>: </p></body></html> - - - false + Variance of selected noise distribution. - - - - - - Inversion time (in ms) for inversion recovery sequences. If 0, no inversion pulse is simulated. + + 10 - 0 + 0.000000000000000 - 999999999 + 999999999.000000000000000 - 1 + 0.001000000000000 - 0 + 50.000000000000000 + + + + + + + Distribution: + + + + + + + Noise distribution + + + Complex Gaussian + + + + + Rician + + + + + + Qt::Horizontal + + + + + + + Qt::Horizontal + + + - + + + true + QFrame::NoFrame QFrame::Raised - + + + 6 + 0 0 0 0 - - 6 - - - + + - <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> + K-Space Line Offset: false - - - - b-value in s/mm² - - - 0 - - - 10000 - - - 100 - - - 1000 - - - - - - - Gradient Directions: - - - - - + + - Number of gradient directions distributed over the half sphere. + A larger offset increases the inensity of the ghost image. - - 0 + + 3 - 10000 + 1.000000000000000 - 1 + 0.010000000000000 - 30 + 0.250000000000000 - - + + + + Add Motion Artifacts + + + false + + + + + + + Add N/2 Ghosts + + + false + + + + + + + Qt::Horizontal + + + + + + + Add ringing artifacts occuring at strong edges in the image. + + + Add Gibbs Ringing + + + false + + + + + + + Qt::Horizontal + + + + + + + Add Noise + + + false + + + + + + + Qt::Horizontal + + + + + + + Add Aliasing + + + false + + + + + + + If > 0, ringing is simulated by by setting the defined percentage of higher frequencies to 0 in k-space. Otherwise, the input to the k-space simulation is generated with twice the resolution and cropped during k-space simulation (much slower). + + + 100 + + + 10 + + + + + + + + + + Qt::Vertical + + + + 20 + 40 + + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + 0 + + + 0 + + + 0 + + + + + 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> + + + + + + Save Parameters + + + + :/QmitkDiffusionImaging/general_icons/download.ico:/QmitkDiffusionImaging/general_icons/download.ico + + + + + + + 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> + + + + - Advanced Options + Load Parameters + + + + :/QmitkDiffusionImaging/general_icons/upload.ico:/QmitkDiffusionImaging/general_icons/upload.ico - - + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + 0 + + + 0 + + + 0 + + + + + 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> + - color: rgb(255, 0, 0); + - Using gradients of selected DWI! + Start Simulation + + + + :/QmitkDiffusionImaging/general_icons/right.ico:/QmitkDiffusionImaging/general_icons/right.ico - - - - QFrame::NoFrame + + + + QGroupBox { + background-color: transparent; +} - - QFrame::Raised + + Input Data - + - 0 + 6 - 0 + 6 - 0 + 6 - 0 + 6 - - - - 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 - - - 20 - - - - - - - Fiber sampling factor which determines the accuracy of the calculated fiber and non-fiber volume fractions. - - - 1 + + + + QFrame::NoFrame - - 1000 + + QFrame::Raised - - 1 + + + 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> - - 20 + + QComboBox::AdjustToMinimumContentsLength - - + + - Fiber sampling factor which determines the accuracy of the calculated fiber and non-fiber volume fractions. + - - 1 + + - - 1000 + + - - 1 + + Fiber Bundle: - - 3 + + false - - - - - - - Use bvals/bvecs files - - - - - - - QFrame::NoFrame - - - QFrame::Raised - - - - 0 - - - 0 - - - 0 - - - 0 - - - - - ... + + + + - - - - - - ... + + + + + - - - - - - + Save path: - + false - + - Bvecs: + 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 + + + + + - Bvals: + 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> - - false + + QComboBox::AdjustToMinimumContentsLength - - - - - - - QGroupBox { - background-color: transparent; -} - - - Inter-axonal Compartment - - - - 6 - - - 6 - - - 6 - - - 6 - - - + + + + true + - Select signal model for intra-axonal compartment. + Stop current simulation. + + + + + + Abort Simulation + + + + :/QmitkDiffusionImaging/general_icons/abort.ico:/QmitkDiffusionImaging/general_icons/abort.ico - - - -- - - - - - Stick Model - - - - - Zeppelin Model - - - - - Tensor Model - - - - - - - - - - - - - - - - - QFrame::NoFrame + + + + Courier + 7 + - - QFrame::Raised + + true - - - 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. - - - - - - - - Qt::Vertical - - - - 20 - 40 - - - - QmitkDataStorageComboBox QComboBox
QmitkDataStorageComboBox.h
 QmitkDataStorageComboBoxWithSelectNone QComboBox
QmitkDataStorageComboBoxWithSelectNone.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 m_FiberBundleComboBox m_MaskComboBox m_TemplateComboBox m_SavePathEdit m_OutputPathButton + m_GenerateImageButton + m_AbortSimulationButton + m_SimulationStatusText + m_LoadParametersButton + m_SaveParametersButton m_SizeX m_SizeY m_SizeZ m_SpacingX m_SpacingY m_SpacingZ + m_UseBvalsBvecsBox + m_LoadBvalsEdit + m_LoadBvalsButton + m_LoadBvecsEdit + m_LoadBvecsButton m_NumGradientsBox m_BvalueBox m_AdvancedOptionsBox_2 + m_AcquisitionTypeBox m_SignalScaleBox + m_NumCoilsBox + m_CoilSensBox m_TEbox m_TRbox + m_TIbox m_LineReadoutTimeBox m_PartialFourier m_T2starBox m_FiberRadius m_ReversePhaseBox m_RelaxationBox m_EnforcePureFiberVoxelsBox m_VolumeFractionsBox m_Compartment1Box m_Comp1VolumeFraction m_Compartment2Box m_Comp2VolumeFraction m_Compartment3Box m_Comp3VolumeFraction m_Compartment4Box m_Comp4VolumeFraction m_AddNoise m_NoiseDistributionBox m_NoiseLevel m_AddSpikes m_SpikeNumBox m_SpikeScaleBox m_AddGhosts m_kOffsetBox m_AddAliasing m_WrapBox m_AddDistortions m_FrequencyMapBox + m_AddDrift + m_DriftFactor m_AddMotion m_RandomMotion + m_MotionVolumesBox m_MaxRotationBoxX m_MaxRotationBoxY m_MaxRotationBoxZ m_MaxTranslationBoxX m_MaxTranslationBoxY m_MaxTranslationBoxZ m_AddEddy m_EddyGradientStrength m_AddGibbsRinging - m_SaveParametersButton - m_LoadParametersButton + m_ZeroRinging