diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkKspaceImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkKspaceImageFilter.cpp index a7916bb498..14d9952c60 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkKspaceImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkKspaceImageFilter.cpp @@ -1,434 +1,451 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkKspaceImageFilter_txx #define __itkKspaceImageFilter_txx //#endif #include #include #include +#include #include "itkKspaceImageFilter.h" #include #include #include #include #include #include namespace itk { template< class ScalarType > KspaceImageFilter< ScalarType >::KspaceImageFilter() : m_Z(0) , m_UseConstantRandSeed(false) , 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]; 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(); outputImage->FillBuffer(m_Spike); m_KSpaceImage = InputImageType::New(); m_KSpaceImage->SetLargestPossibleRegion( region ); m_KSpaceImage->SetBufferedRegion( region ); m_KSpaceImage->SetRequestedRegion( region ); m_KSpaceImage->Allocate(); m_KSpaceImage->FillBuffer(0.0); 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); 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; } } + + 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: + noiseVar = m_Parameters->m_SignalGen.m_PartialFourier*m_Parameters->m_SignalGen.m_NoiseVariance/(kyMax*kxMax); } 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) { itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer randGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); randGen->SetSeed(); if (m_UseConstantRandSeed) // always generate the same random numbers? { randGen->SetSeed(0); } else { randGen->SetSeed(); } typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); typedef ImageRegionConstIterator< InputImageType > InputIteratorType; - float kxMax = m_Parameters->m_SignalGen.m_CroppedRegion.GetSize(0); - float kyMax = m_Parameters->m_SignalGen.m_CroppedRegion.GetSize(1); - float xMax = m_CompartmentImages.at(0)->GetLargestPossibleRegion().GetSize(0); // scanner coverage in x-direction - float yMax = m_CompartmentImages.at(0)->GetLargestPossibleRegion().GetSize(1); // scanner coverage in y-direction - float 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) - - float numPix = kxMax*kyMax; - // 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); - 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) + { + ++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) + { + vcl_complex zero = vcl_complex(0, 0); + oit.Set(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) * (1.0-std::exp(-(m_Parameters->m_SignalGen.m_tRep + tRf)/m_T1[i])) ); } - // get current k-space index (depends on the chosen k-space readout scheme) - itk::Index< 2 > kIdx = m_ReadoutScheme->GetActualKspaceIndex(out_idx); + // 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; + } - // partial fourier - bool partial_fourier_skip = false; - if (kIdx[1]>kyMax*m_Parameters->m_SignalGen.m_PartialFourier) - partial_fourier_skip = true; + // pull stuff out of inner loop + t /= 1000; + kx /= xMax; + ky /= yMaxFov; - if (!partial_fourier_skip) + // 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() ) { - // 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; + typename InputImageType::IndexType input_idx = it.GetIndex(); - if (static_cast(kyMax)%2==1) - ky -= (kyMax-1)/2; + // 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 - ky -= kyMax/2; + x -= xMax/2; + if (static_cast(yMax)%2==1) + y -= (yMax-1)/2; + else + y -= yMax/2; - // 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; + // 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 - kx += m_Parameters->m_SignalGen.m_KspaceLineOffset; - } + f_real += m_CompartmentImages[i]->GetPixel(input_idx) * m_Parameters->m_SignalGen.m_SignalScale; - // pull stuff out of inner loop - t /= 1000; - kx /= xMax; - ky /= yMaxFov; - auto yMaxFov_half = yMaxFov/2; - - // 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() ) + // 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) { - 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; - } + 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); + 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 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) + // 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()) { - 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); - } + 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; - } + // 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 )) ); + // 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; + ++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; + 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; - if (m_Parameters->m_SignalGen.m_NoiseVariance>0 && m_Parameters->m_Misc.m_DoAddNoise) - s = vcl_complex(s.real()+randGen->GetNormalVariate(0,noiseVar), s.imag()+randGen->GetNormalVariate(0,noiseVar)); + if (m_Parameters->m_SignalGen.m_NoiseVariance>0 && m_Parameters->m_Misc.m_DoAddNoise) + s = vcl_complex(s.real()+randGen->GetNormalVariate(0,noiseVar), s.imag()+randGen->GetNormalVariate(0,noiseVar)); - outputImage->SetPixel(kIdx, s); - m_KSpaceImage->SetPixel(kIdx, sqrt(s.imag()*s.imag()+s.real()*s.real()) ); - } + 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; } itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer randGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); randGen->SetSeed(); if (m_UseConstantRandSeed) // always generate the same random numbers? randGen->SetSeed(0); else randGen->SetSeed(); m_Spike *= m_Parameters->m_SignalGen.m_SpikeAmplitude; itk::Index< 2 > spikeIdx; for (unsigned int i=0; iGetIntegerVariate()%kxMax; spikeIdx[1] = 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/itkKspaceImageFilter.h b/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkKspaceImageFilter.h index 519b8f4234..69649e6af5 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkKspaceImageFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkKspaceImageFilter.h @@ -1,146 +1,156 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ /*=================================================================== This file is based heavily on a corresponding ITK filter. ===================================================================*/ #ifndef __itkKspaceImageFilter_h_ #define __itkKspaceImageFilter_h_ #include #include #include #include #include #include #include #include namespace itk{ /** * \brief Simulates k-space acquisition of one slice with a single shot EPI sequence. Enables the simulation of various effects occuring during real MR acquisitions: * - T2 signal relaxation * - Spikes * - N/2 Ghosts * - Aliasing (wrap around) * - Image distortions (off-frequency effects) * - Gibbs ringing * - Eddy current effects * Based on a discrete fourier transformation. * See "Fiberfox: Facilitating the creation of realistic white matter software phantoms" (DOI: 10.1002/mrm.25045) for details. */ template< class ScalarType > class KspaceImageFilter : public ImageSource< Image< vcl_complex< ScalarType >, 2 > > { public: typedef KspaceImageFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageSource< Image< vcl_complex< ScalarType >, 2 > > Superclass; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(KspaceImageFilter, ImageToImageFilter) typedef typename itk::Image< ScalarType, 2 > InputImageType; typedef typename InputImageType::Pointer InputImagePointerType; typedef typename Superclass::OutputImageType OutputImageType; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; typedef itk::Matrix MatrixType; typedef itk::Point Point2D; typedef itk::Vector< float,3> VectorType; itkSetMacro( SpikesPerSlice, unsigned int ) ///< Number of spikes per slice. Corresponding parameter in fiberfox parameter object specifies the number of spikes for the whole image and can thus not be used here. itkSetMacro( Z, double ) ///< Slice position, necessary for eddy current simulation. itkSetMacro( UseConstantRandSeed, bool ) ///< Use constant seed for random generator for reproducible results. ONLY USE FOR TESTING PURPOSES! itkSetMacro( Translation, VectorType ) itkSetMacro( RotationMatrix, MatrixType ) itkSetMacro( Zidx, int ) itkSetMacro( FiberBundle, FiberBundle::Pointer ) itkSetMacro( CoilPosition, VectorType ) itkGetMacro( KSpaceImage, typename InputImageType::Pointer ) ///< k-space magnitude image itkGetMacro( SpikeLog, std::string ) void SetParameters( FiberfoxParameters* param ){ m_Parameters = param; } void SetCompartmentImages( std::vector< InputImagePointerType > cImgs ) { m_CompartmentImages=cImgs; } ///< One signal image per compartment. void SetT2( std::vector< float > t2Vector ) { m_T2=t2Vector; } ///< One T2 relaxation constant per compartment image. void SetT1( std::vector< float > t1Vector ) { m_T1=t1Vector; } ///< One T1 relaxation constant per compartment image. void SetDiffusionGradientDirection(itk::Vector g) { m_DiffusionGradientDirection=g; } ///< Gradient direction is needed for eddy current simulation. protected: KspaceImageFilter(); ~KspaceImageFilter() override {} float CoilSensitivity(VectorType& pos); void BeforeThreadedGenerateData() override; void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType threadID) override; void AfterThreadedGenerateData() override; VectorType m_CoilPosition; FiberfoxParameters* m_Parameters; std::vector< float > m_T2; std::vector< float > m_T1; std::vector< InputImagePointerType > m_CompartmentImages; itk::Vector m_DiffusionGradientDirection; float m_Z; int m_Zidx; bool m_UseConstantRandSeed; unsigned int m_SpikesPerSlice; FiberBundle::Pointer m_FiberBundle; float m_Gamma; VectorType m_Translation; ///< used to find correct point in frequency map (head motion) MatrixType m_RotationMatrix; float m_TransX; float m_TransY; float m_TransZ; bool m_IsBaseline; vcl_complex m_Spike; MatrixType m_Transform; std::string m_SpikeLog; float m_CoilSensitivityFactor; typename InputImageType::Pointer m_KSpaceImage; typename InputImageType::Pointer m_TimeFromEchoImage; typename InputImageType::Pointer m_ReadoutTimeImage; AcquisitionType* m_ReadoutScheme; typename itk::Image< ScalarType, 2 >::Pointer m_MovedFmap; + int ringing_lines_x; + int ringing_lines_y; + float kxMax; + float kyMax; + float xMax; + float yMax; + float yMaxFov; + float yMaxFov_half; + float numPix; + float noiseVar; private: }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkKspaceImageFilter.cpp" #endif #endif //__itkKspaceImageFilter_h_ diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp index 24949e8598..4111beb43c 100755 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp @@ -1,1770 +1,1775 @@ /*=================================================================== 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 namespace itk { template< class PixelType > TractsToDWIImageFilter< PixelType >::TractsToDWIImageFilter() : m_StatusText("") , m_UseConstantRandSeed(false) , m_RandGen(itk::Statistics::MersenneTwisterRandomVariateGenerator::New()) { m_RandGen->SetSeed(); m_DoubleInterpolator = itk::LinearInterpolateImageFunction< ItkDoubleImgType, float >::New(); m_NullDir.Fill(0); } template< class PixelType > TractsToDWIImageFilter< PixelType >::~TractsToDWIImageFilter() { } template< class PixelType > TractsToDWIImageFilter< PixelType >::DoubleDwiType::Pointer TractsToDWIImageFilter< PixelType >:: SimulateKspaceAcquisition( std::vector< DoubleDwiType::Pointer >& compartment_images ) { unsigned int numFiberCompartments = m_Parameters.m_FiberModelList.size(); // create slice object ImageRegion<2> sliceRegion; sliceRegion.SetSize(0, m_WorkingImageRegion.GetSize()[0]); sliceRegion.SetSize(1, m_WorkingImageRegion.GetSize()[1]); Vector< double, 2 > sliceSpacing; sliceSpacing[0] = m_WorkingSpacing[0]; sliceSpacing[1] = m_WorkingSpacing[1]; DoubleDwiType::PixelType nullPix; nullPix.SetSize(compartment_images.at(0)->GetVectorLength()); nullPix.Fill(0.0); auto magnitudeDwiImage = DoubleDwiType::New(); magnitudeDwiImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); magnitudeDwiImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin ); magnitudeDwiImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); magnitudeDwiImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); magnitudeDwiImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); magnitudeDwiImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); magnitudeDwiImage->SetVectorLength( compartment_images.at(0)->GetVectorLength() ); magnitudeDwiImage->Allocate(); magnitudeDwiImage->FillBuffer(nullPix); m_PhaseImage = DoubleDwiType::New(); m_PhaseImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); m_PhaseImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin ); m_PhaseImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); m_PhaseImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_PhaseImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_PhaseImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_PhaseImage->SetVectorLength( compartment_images.at(0)->GetVectorLength() ); m_PhaseImage->Allocate(); m_PhaseImage->FillBuffer(nullPix); m_KspaceImage = DoubleDwiType::New(); m_KspaceImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); m_KspaceImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin ); m_KspaceImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); m_KspaceImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_KspaceImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_KspaceImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_KspaceImage->SetVectorLength( m_Parameters.m_SignalGen.m_NumberOfCoils ); m_KspaceImage->Allocate(); m_KspaceImage->FillBuffer(nullPix); std::list< unsigned int > spikeVolume; if (m_Parameters.m_Misc.m_DoAddSpikes) for (unsigned int i=0; iGetIntegerVariate()%(compartment_images.at(0)->GetVectorLength())); // 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 (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 max_threads = omp_get_max_threads(); int out_threads = Math::ceil(std::sqrt(max_threads)); int in_threads = Math::floor(std::sqrt(max_threads)); PrintToLog("Parallel volumes: " + boost::lexical_cast(out_threads), false, true, true); PrintToLog("Threads per slice: " + boost::lexical_cast(in_threads), false, true, true); PrintToLog("0% 10 20 30 40 50 60 70 80 90 100%", false, true, false); PrintToLog("|----|----|----|----|----|----|----|----|----|----|\n*", false, false, false); unsigned long lastTick = 0; boost::progress_display disp(compartment_images.at(0)->GetVectorLength()*compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2)); #pragma omp parallel for num_threads(out_threads) for (int g=0; g(compartment_images.at(0)->GetVectorLength()); g++) { if (this->GetAbortGenerateData()) continue; std::list< unsigned int > spikeSlice; #pragma omp critical { for (auto sv : spikeVolume) if (sv == static_cast(g)) spikeSlice.push_back(m_RandGen->GetIntegerVariate()%compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2)); spikeVolume.remove(g); } for (unsigned int z=0; zGetLargestPossibleRegion().GetSize(2); z++) { std::vector< Float2DImageType::Pointer > 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()); } int numSpikes = 0; for (auto ss : spikeSlice) if (ss == z) ++numSpikes; spikeSlice.remove(z); int spikeCoil = m_RandGen->GetIntegerVariate()%m_Parameters.m_SignalGen.m_NumberOfCoils; if (this->GetAbortGenerateData()) continue; for (int c=0; c::New(); idft->SetCompartmentImages(compartment_slices); idft->SetT2(t2Vector); idft->SetT1(t1Vector); idft->SetUseConstantRandSeed(m_UseConstantRandSeed); idft->SetParameters(&m_Parameters); idft->SetZ((float)z-(float)( compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2) -compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2)%2 ) / 2.0); idft->SetZidx(z); idft->SetCoilPosition(coilPositions.at(c)); idft->SetFiberBundle(m_FiberBundle); idft->SetTranslation(m_Translations.at(g)); idft->SetRotationMatrix(m_RotationsInv.at(g)); idft->SetDiffusionGradientDirection(m_Parameters.m_SignalGen.GetGradientDirection(g)); if (c==spikeCoil) idft->SetSpikesPerSlice(numSpikes); idft->SetNumberOfThreads(in_threads); idft->Update(); #pragma omp critical if (c==spikeCoil && 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 (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) + 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) + 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: { 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_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) - PrintToLog("Simulating ringing artifacts", false); + { + 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.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.cpp index ba8d741ac7..8292c1f681 100644 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.cpp +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.cpp @@ -1,1077 +1,1079 @@ /*=================================================================== 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. ===================================================================*/ #define RAPIDXML_NO_EXCEPTIONS #include #include #include #include #include #include #include #include #include mitk::FiberfoxParameters::FiberfoxParameters() : m_NoiseModel(nullptr) { mitk::StickModel* model_aniso = new mitk::StickModel(); model_aniso->m_CompartmentId = 1; m_FiberModelList.push_back(model_aniso); mitk::BallModel* model_iso = new mitk::BallModel(); model_iso->m_CompartmentId = 3; m_NonFiberModelList.push_back(model_iso); } mitk::FiberfoxParameters::FiberfoxParameters(const mitk::FiberfoxParameters& params) : m_NoiseModel(nullptr) { m_FiberGen = params.m_FiberGen; m_SignalGen = params.m_SignalGen; m_Misc = params.m_Misc; if (params.m_NoiseModel!=nullptr) { if (dynamic_cast*>(params.m_NoiseModel.get())) m_NoiseModel = std::make_shared< mitk::RicianNoiseModel<> >(); else if (dynamic_cast*>(params.m_NoiseModel.get())) m_NoiseModel = std::make_shared< mitk::ChiSquareNoiseModel<> >(); m_NoiseModel->SetNoiseVariance(params.m_NoiseModel->GetNoiseVariance()); } for (unsigned int i=0; i* outModel = nullptr; mitk::DiffusionSignalModel<>* signalModel = nullptr; if (i*>(signalModel)) outModel = new mitk::StickModel<>(dynamic_cast*>(signalModel)); else if (dynamic_cast*>(signalModel)) outModel = new mitk::TensorModel<>(dynamic_cast*>(signalModel)); else if (dynamic_cast*>(signalModel)) outModel = new mitk::RawShModel<>(dynamic_cast*>(signalModel)); else if (dynamic_cast*>(signalModel)) outModel = new mitk::BallModel<>(dynamic_cast*>(signalModel)); else if (dynamic_cast*>(signalModel)) outModel = new mitk::AstroStickModel<>(dynamic_cast*>(signalModel)); else if (dynamic_cast*>(signalModel)) outModel = new mitk::DotModel<>(dynamic_cast*>(signalModel)); if (i theta; theta.set_size(NPoints); vnl_vector phi; phi.set_size(NPoints); double C = sqrt(4*itk::Math::pi); phi(0) = 0.0; phi(NPoints-1) = 0.0; for(int i=0; i0 && i mitk::SignalGenerationParameters::GetBaselineIndices() { std::vector< int > result; for( unsigned int i=0; im_GradientDirections.size(); i++) if (m_GradientDirections.at(i).GetNorm()<0.0001) result.push_back(i); return result; } unsigned int mitk::SignalGenerationParameters::GetFirstBaselineIndex() { for( unsigned int i=0; im_GradientDirections.size(); i++) if (m_GradientDirections.at(i).GetNorm()<0.0001) return i; return -1; } bool mitk::SignalGenerationParameters::IsBaselineIndex(unsigned int idx) { if (m_GradientDirections.size()>idx && m_GradientDirections.at(idx).GetNorm()<0.0001) return true; return false; } unsigned int mitk::SignalGenerationParameters::GetNumWeightedVolumes() { return m_NumGradients; } unsigned int mitk::SignalGenerationParameters::GetNumBaselineVolumes() { return m_NumBaseline; } unsigned int mitk::SignalGenerationParameters::GetNumVolumes() { return m_GradientDirections.size(); } mitk::SignalGenerationParameters::GradientListType mitk::SignalGenerationParameters::GetGradientDirections() { return m_GradientDirections; } mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer mitk::SignalGenerationParameters::GetItkGradientContainer() { int c = 0; mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer out = mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::New(); for (auto g : m_GradientDirections) { mitk::DiffusionPropertyHelper::GradientDirectionType vnl_dir; vnl_dir[0] = g[0]; vnl_dir[1] = g[1]; vnl_dir[2] = g[2]; out->InsertElement(c, vnl_dir); ++c; } return out; } mitk::SignalGenerationParameters::GradientType mitk::SignalGenerationParameters::GetGradientDirection(unsigned int i) { return m_GradientDirections.at(i); } void mitk::SignalGenerationParameters::SetNumWeightedVolumes(int numGradients) { m_NumGradients = numGradients; GenerateGradientHalfShell(); } std::vector< int > mitk::SignalGenerationParameters::GetBvalues() { std::vector< int > bVals; for( GradientType g : m_GradientDirections) { float norm = g.GetNorm(); int bVal = std::round(norm*norm*m_Bvalue); if ( std::find(bVals.begin(), bVals.end(), bVal) == bVals.end() ) bVals.push_back(bVal); } return bVals; } double mitk::SignalGenerationParameters::GetBvalue() { return m_Bvalue; } void mitk::SignalGenerationParameters::SetGradienDirections(GradientListType gradientList) { m_GradientDirections = gradientList; m_NumGradients = 0; m_NumBaseline = 0; for( unsigned int i=0; im_GradientDirections.size(); i++) { float norm = m_GradientDirections.at(i).GetNorm(); if (norm>0.0001) m_NumGradients++; else m_NumBaseline++; } } void mitk::SignalGenerationParameters::SetGradienDirections(mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradientList) { m_NumGradients = 0; m_NumBaseline = 0; m_GradientDirections.clear(); for( unsigned int i=0; iSize(); i++) { GradientType g; g[0] = gradientList->at(i)[0]; g[1] = gradientList->at(i)[1]; g[2] = gradientList->at(i)[2]; m_GradientDirections.push_back(g); float norm = m_GradientDirections.at(i).GetNorm(); if (norm>0.0001) m_NumGradients++; else m_NumBaseline++; } } void mitk::SignalGenerationParameters::ApplyDirectionMatrix() { auto imageRotationMatrix = m_ImageDirection.GetVnlMatrix(); GradientListType rotated_gradients; for(auto g : m_GradientDirections) { vnl_vector vec = g.GetVnlVector(); vec = vec.pre_multiply(imageRotationMatrix); GradientType g2; g2[0] = vec[0]; g2[1] = vec[1]; g2[2] = vec[2]; rotated_gradients.push_back(g2); } m_GradientDirections = rotated_gradients; } void mitk::FiberfoxParameters::ApplyDirectionMatrix() { m_SignalGen.ApplyDirectionMatrix(); UpdateSignalModels(); } void mitk::FiberfoxParameters::SaveParameters(std::string filename) { if(filename.empty()) return; if(".ffp"!=filename.substr(filename.size()-4, 4)) filename += ".ffp"; const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, nullptr ); if ( locale.compare(currLocale)!=0 ) { try { setlocale(LC_ALL, locale.c_str()); } catch(...) { MITK_INFO << "Could not set locale " << locale; } } boost::property_tree::ptree parameters; // fiber generation parameters parameters.put("fiberfox.fibers.distribution", m_FiberGen.m_Distribution); parameters.put("fiberfox.fibers.variance", m_FiberGen.m_Variance); parameters.put("fiberfox.fibers.density", m_FiberGen.m_Density); parameters.put("fiberfox.fibers.spline.sampling", m_FiberGen.m_Sampling); parameters.put("fiberfox.fibers.spline.tension", m_FiberGen.m_Tension); parameters.put("fiberfox.fibers.spline.continuity", m_FiberGen.m_Continuity); parameters.put("fiberfox.fibers.spline.bias", m_FiberGen.m_Bias); parameters.put("fiberfox.fibers.rotation.x", m_FiberGen.m_Rotation[0]); parameters.put("fiberfox.fibers.rotation.y", m_FiberGen.m_Rotation[1]); parameters.put("fiberfox.fibers.rotation.z", m_FiberGen.m_Rotation[2]); parameters.put("fiberfox.fibers.translation.x", m_FiberGen.m_Translation[0]); parameters.put("fiberfox.fibers.translation.y", m_FiberGen.m_Translation[1]); parameters.put("fiberfox.fibers.translation.z", m_FiberGen.m_Translation[2]); parameters.put("fiberfox.fibers.scale.x", m_FiberGen.m_Scale[0]); parameters.put("fiberfox.fibers.scale.y", m_FiberGen.m_Scale[1]); parameters.put("fiberfox.fibers.scale.z", m_FiberGen.m_Scale[2]); // image generation parameters parameters.put("fiberfox.image.basic.size.x", m_SignalGen.m_ImageRegion.GetSize(0)); parameters.put("fiberfox.image.basic.size.y", m_SignalGen.m_ImageRegion.GetSize(1)); parameters.put("fiberfox.image.basic.size.z", m_SignalGen.m_ImageRegion.GetSize(2)); parameters.put("fiberfox.image.basic.spacing.x", m_SignalGen.m_ImageSpacing[0]); parameters.put("fiberfox.image.basic.spacing.y", m_SignalGen.m_ImageSpacing[1]); parameters.put("fiberfox.image.basic.spacing.z", m_SignalGen.m_ImageSpacing[2]); parameters.put("fiberfox.image.basic.origin.x", m_SignalGen.m_ImageOrigin[0]); parameters.put("fiberfox.image.basic.origin.y", m_SignalGen.m_ImageOrigin[1]); parameters.put("fiberfox.image.basic.origin.z", m_SignalGen.m_ImageOrigin[2]); parameters.put("fiberfox.image.basic.direction.d1", m_SignalGen.m_ImageDirection[0][0]); parameters.put("fiberfox.image.basic.direction.d2", m_SignalGen.m_ImageDirection[0][1]); parameters.put("fiberfox.image.basic.direction.d3", m_SignalGen.m_ImageDirection[0][2]); parameters.put("fiberfox.image.basic.direction.d4", m_SignalGen.m_ImageDirection[1][0]); parameters.put("fiberfox.image.basic.direction.d5", m_SignalGen.m_ImageDirection[1][1]); parameters.put("fiberfox.image.basic.direction.d6", m_SignalGen.m_ImageDirection[1][2]); parameters.put("fiberfox.image.basic.direction.d7", m_SignalGen.m_ImageDirection[2][0]); parameters.put("fiberfox.image.basic.direction.d8", m_SignalGen.m_ImageDirection[2][1]); parameters.put("fiberfox.image.basic.direction.d9", m_SignalGen.m_ImageDirection[2][2]); mitk::gradients::WriteBvalsBvecs(filename+".bvals", filename+".bvecs", m_SignalGen.GetItkGradientContainer(), m_SignalGen.m_Bvalue); parameters.put("fiberfox.image.acquisitiontype", m_SignalGen.m_AcquisitionType); parameters.put("fiberfox.image.coilsensitivityprofile", m_SignalGen.m_CoilSensitivityProfile); parameters.put("fiberfox.image.numberofcoils", m_SignalGen.m_NumberOfCoils); parameters.put("fiberfox.image.reversephase", m_SignalGen.m_ReversePhase); parameters.put("fiberfox.image.partialfourier", m_SignalGen.m_PartialFourier); parameters.put("fiberfox.image.noisevariance", m_SignalGen.m_NoiseVariance); parameters.put("fiberfox.image.trep", m_SignalGen.m_tRep); parameters.put("fiberfox.image.signalScale", m_SignalGen.m_SignalScale); parameters.put("fiberfox.image.tEcho", m_SignalGen.m_tEcho); parameters.put("fiberfox.image.tLine", m_SignalGen.m_tLine); parameters.put("fiberfox.image.tInhom", m_SignalGen.m_tInhom); parameters.put("fiberfox.image.simulatekspace", m_SignalGen.m_SimulateKspaceAcquisition); parameters.put("fiberfox.image.axonRadius", m_SignalGen.m_AxonRadius); parameters.put("fiberfox.image.doSimulateRelaxation", m_SignalGen.m_DoSimulateRelaxation); parameters.put("fiberfox.image.doDisablePartialVolume", m_SignalGen.m_DoDisablePartialVolume); parameters.put("fiberfox.image.artifacts.spikesnum", m_SignalGen.m_Spikes); parameters.put("fiberfox.image.artifacts.spikesscale", m_SignalGen.m_SpikeAmplitude); parameters.put("fiberfox.image.artifacts.kspaceLineOffset", m_SignalGen.m_KspaceLineOffset); parameters.put("fiberfox.image.artifacts.eddyStrength", m_SignalGen.m_EddyStrength); parameters.put("fiberfox.image.artifacts.eddyTau", m_SignalGen.m_Tau); parameters.put("fiberfox.image.artifacts.aliasingfactor", m_SignalGen.m_CroppingFactor); parameters.put("fiberfox.image.artifacts.drift", m_SignalGen.m_Drift); parameters.put("fiberfox.image.artifacts.doAddMotion", m_SignalGen.m_DoAddMotion); parameters.put("fiberfox.image.artifacts.randomMotion", m_SignalGen.m_DoRandomizeMotion); parameters.put("fiberfox.image.artifacts.translation0", m_SignalGen.m_Translation[0]); parameters.put("fiberfox.image.artifacts.translation1", m_SignalGen.m_Translation[1]); parameters.put("fiberfox.image.artifacts.translation2", m_SignalGen.m_Translation[2]); parameters.put("fiberfox.image.artifacts.rotation0", m_SignalGen.m_Rotation[0]); parameters.put("fiberfox.image.artifacts.rotation1", m_SignalGen.m_Rotation[1]); parameters.put("fiberfox.image.artifacts.rotation2", m_SignalGen.m_Rotation[2]); parameters.put("fiberfox.image.artifacts.motionvolumes", m_Misc.m_MotionVolumesBox); parameters.put("fiberfox.image.artifacts.addringing", m_SignalGen.m_DoAddGibbsRinging); + parameters.put("fiberfox.image.artifacts.zeroringing", m_SignalGen.m_ZeroRinging); parameters.put("fiberfox.image.artifacts.addnoise", m_Misc.m_DoAddNoise); parameters.put("fiberfox.image.artifacts.addghosts", m_Misc.m_DoAddGhosts); parameters.put("fiberfox.image.artifacts.addaliasing", m_Misc.m_DoAddAliasing); parameters.put("fiberfox.image.artifacts.addspikes", m_Misc.m_DoAddSpikes); parameters.put("fiberfox.image.artifacts.addeddycurrents", m_Misc.m_DoAddEddyCurrents); parameters.put("fiberfox.image.artifacts.doAddDistortions", m_Misc.m_DoAddDistortions); parameters.put("fiberfox.image.artifacts.doAddDrift", m_SignalGen.m_DoAddDrift); parameters.put("fiberfox.image.outputvolumefractions", m_Misc.m_CheckOutputVolumeFractionsBox); parameters.put("fiberfox.image.showadvanced", m_Misc.m_CheckAdvancedSignalOptionsBox); parameters.put("fiberfox.image.signalmodelstring", m_Misc.m_SignalModelString); parameters.put("fiberfox.image.artifactmodelstring", m_Misc.m_ArtifactModelString); parameters.put("fiberfox.image.outpath", m_Misc.m_OutputPath); parameters.put("fiberfox.fibers.realtime", m_Misc.m_CheckRealTimeFibersBox); parameters.put("fiberfox.fibers.showadvanced", m_Misc.m_CheckAdvancedFiberOptionsBox); parameters.put("fiberfox.fibers.constantradius", m_Misc.m_CheckConstantRadiusBox); parameters.put("fiberfox.fibers.includeFiducials", m_Misc.m_CheckIncludeFiducialsBox); if (m_NoiseModel!=nullptr) { parameters.put("fiberfox.image.artifacts.noisevariance", m_NoiseModel->GetNoiseVariance()); if (dynamic_cast*>(m_NoiseModel.get())) parameters.put("fiberfox.image.artifacts.noisetype", "rice"); else if (dynamic_cast*>(m_NoiseModel.get())) parameters.put("fiberfox.image.artifacts.noisetype", "chisquare"); } for (std::size_t i=0; i* signalModel = nullptr; if (i(i)+".type", "fiber"); } else { signalModel = m_NonFiberModelList.at(i-m_FiberModelList.size()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".type", "non-fiber"); } if (dynamic_cast*>(signalModel)) { mitk::StickModel<>* model = dynamic_cast*>(signalModel); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".model", "stick"); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".d", model->GetDiffusivity()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".t2", model->GetT2()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".t1", model->GetT1()); } else if (dynamic_cast*>(signalModel)) { mitk::TensorModel<>* model = dynamic_cast*>(signalModel); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".model", "tensor"); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".d1", model->GetDiffusivity1()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".d2", model->GetDiffusivity2()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".d3", model->GetDiffusivity3()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".t2", model->GetT2()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".t1", model->GetT1()); } else if (dynamic_cast*>(signalModel)) { mitk::RawShModel<>* model = dynamic_cast*>(signalModel); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".model", "prototype"); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".minFA", model->GetFaRange().first); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".maxFA", model->GetFaRange().second); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".minADC", model->GetAdcRange().first); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".maxADC", model->GetAdcRange().second); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".maxNumSamples", model->GetMaxNumKernels()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".numSamples", model->GetNumberOfKernels()); int shOrder = model->GetShOrder(); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".numCoeffs", (shOrder*shOrder + shOrder + 2)/2 + shOrder); for (unsigned int j=0; jGetNumberOfKernels(); j++) { vnl_vector< double > coeffs = model->GetCoefficients(j); for (unsigned int k=0; k(i)+".kernels."+boost::lexical_cast(j)+".coeffs."+boost::lexical_cast(k), coeffs[k]); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".kernels."+boost::lexical_cast(j)+".B0", model->GetBaselineSignal(j)); } } else if (dynamic_cast*>(signalModel)) { mitk::BallModel<>* model = dynamic_cast*>(signalModel); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".model", "ball"); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".d", model->GetDiffusivity()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".t2", model->GetT2()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".t1", model->GetT1()); } else if (dynamic_cast*>(signalModel)) { mitk::AstroStickModel<>* model = dynamic_cast*>(signalModel); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".model", "astrosticks"); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".d", model->GetDiffusivity()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".t2", model->GetT2()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".t1", model->GetT1()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".randomize", model->GetRandomizeSticks()); } else if (dynamic_cast*>(signalModel)) { mitk::DotModel<>* model = dynamic_cast*>(signalModel); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".model", "dot"); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".t2", model->GetT2()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".t1", model->GetT1()); } if (signalModel!=nullptr) { parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".ID", signalModel->m_CompartmentId); if (signalModel->GetVolumeFractionImage().IsNotNull()) { try{ itk::ImageFileWriter::Pointer writer = itk::ImageFileWriter::New(); writer->SetFileName(filename+"_VOLUME"+boost::lexical_cast(signalModel->m_CompartmentId)+".nii.gz"); writer->SetInput(signalModel->GetVolumeFractionImage()); writer->Update(); MITK_INFO << "Volume fraction image for compartment "+boost::lexical_cast(signalModel->m_CompartmentId)+" saved."; } catch(...) { } } } } boost::property_tree::xml_writer_settings writerSettings(' ', 2); boost::property_tree::xml_parser::write_xml(filename, parameters, std::locale(), writerSettings); try{ itk::ImageFileWriter::Pointer writer = itk::ImageFileWriter::New(); writer->SetFileName(filename+"_FMAP.nii.gz"); writer->SetInput(m_SignalGen.m_FrequencyMap); writer->Update(); } catch(...) { MITK_INFO << "No frequency map saved."; } try{ itk::ImageFileWriter::Pointer writer = itk::ImageFileWriter::New(); writer->SetFileName(filename+"_MASK.nii.gz"); writer->SetInput(m_SignalGen.m_MaskImage); writer->Update(); } catch(...) { MITK_INFO << "No mask image saved."; } setlocale(LC_ALL, currLocale.c_str()); } template< class ParameterType > ParameterType mitk::FiberfoxParameters::ReadVal(boost::property_tree::ptree::value_type const& v, std::string tag, ParameterType defaultValue, bool essential) { try { return v.second.get(tag); } catch (...) { if (essential) { mitkThrow() << "Parameter file corrupted. Essential tag is missing: '" << tag << "'"; } if (tag!="artifacts.noisetype") { MITK_INFO << "Tag '" << tag << "' not found. Using default value '" << defaultValue << "'."; m_MissingTags += "\n- "; m_MissingTags += tag; } return defaultValue; } } void mitk::FiberfoxParameters::UpdateSignalModels() { for (mitk::DiffusionSignalModel<>* m : m_FiberModelList) { m->SetGradientList(m_SignalGen.m_GradientDirections); m->SetBvalue(m_SignalGen.m_Bvalue); } for (mitk::DiffusionSignalModel<>* m : m_NonFiberModelList) { m->SetGradientList(m_SignalGen.m_GradientDirections); m->SetBvalue(m_SignalGen.m_Bvalue); } } void mitk::FiberfoxParameters::SetNumWeightedVolumes(int numGradients) { m_SignalGen.SetNumWeightedVolumes(numGradients); UpdateSignalModels(); } void mitk::FiberfoxParameters::SetGradienDirections(mitk::SignalGenerationParameters::GradientListType gradientList) { m_SignalGen.SetGradienDirections(gradientList); UpdateSignalModels(); } void mitk::FiberfoxParameters::SetGradienDirections(mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradientList) { m_SignalGen.SetGradienDirections(gradientList); UpdateSignalModels(); } void mitk::FiberfoxParameters::SetBvalue(double Bvalue) { m_SignalGen.m_Bvalue = Bvalue; UpdateSignalModels(); } void mitk::FiberfoxParameters::GenerateGradientHalfShell() { m_SignalGen.GenerateGradientHalfShell(); UpdateSignalModels(); } void mitk::FiberfoxParameters::LoadParameters(std::string filename, bool fix_seed) { itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer randgen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); if (fix_seed) { srand(0); randgen->SetSeed(0); } else { srand(time(0)); randgen->SetSeed(); } m_MissingTags = ""; if(filename.empty()) { return; } const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, nullptr ); if ( locale.compare(currLocale)!=0 ) { try { setlocale(LC_ALL, locale.c_str()); } catch(...) { MITK_INFO << "Could not set locale " << locale; } } boost::property_tree::ptree parameterTree; boost::property_tree::xml_parser::read_xml( filename, parameterTree ); m_FiberModelList.clear(); m_NonFiberModelList.clear(); if (m_NoiseModel) { m_NoiseModel = nullptr; } BOOST_FOREACH( boost::property_tree::ptree::value_type const& v1, parameterTree.get_child("fiberfox") ) { if( v1.first == "fibers" ) { m_Misc.m_CheckRealTimeFibersBox = ReadVal(v1,"realtime", m_Misc.m_CheckRealTimeFibersBox); m_Misc.m_CheckAdvancedFiberOptionsBox = ReadVal(v1,"showadvanced", m_Misc.m_CheckAdvancedFiberOptionsBox); m_Misc.m_CheckConstantRadiusBox = ReadVal(v1,"constantradius", m_Misc.m_CheckConstantRadiusBox); m_Misc.m_CheckIncludeFiducialsBox = ReadVal(v1,"includeFiducials", m_Misc.m_CheckIncludeFiducialsBox); switch (ReadVal(v1,"distribution", 0)) { case 0: m_FiberGen.m_Distribution = FiberGenerationParameters::DISTRIBUTE_UNIFORM; break; case 1: m_FiberGen.m_Distribution = FiberGenerationParameters::DISTRIBUTE_GAUSSIAN; break; default: m_FiberGen.m_Distribution = FiberGenerationParameters::DISTRIBUTE_UNIFORM; } m_FiberGen.m_Variance = ReadVal(v1,"variance", m_FiberGen.m_Variance); m_FiberGen.m_Density = ReadVal(v1,"density", m_FiberGen.m_Density); m_FiberGen.m_Sampling = ReadVal(v1,"spline.sampling", m_FiberGen.m_Sampling); m_FiberGen.m_Tension = ReadVal(v1,"spline.tension", m_FiberGen.m_Tension); m_FiberGen.m_Continuity = ReadVal(v1,"spline.continuity", m_FiberGen.m_Continuity); m_FiberGen.m_Bias = ReadVal(v1,"spline.bias", m_FiberGen.m_Bias); m_FiberGen.m_Rotation[0] = ReadVal(v1,"rotation.x", m_FiberGen.m_Rotation[0]); m_FiberGen.m_Rotation[1] = ReadVal(v1,"rotation.y", m_FiberGen.m_Rotation[1]); m_FiberGen.m_Rotation[2] = ReadVal(v1,"rotation.z", m_FiberGen.m_Rotation[2]); m_FiberGen.m_Translation[0] = ReadVal(v1,"translation.x", m_FiberGen.m_Translation[0]); m_FiberGen.m_Translation[1] = ReadVal(v1,"translation.y", m_FiberGen.m_Translation[1]); m_FiberGen.m_Translation[2] = ReadVal(v1,"translation.z", m_FiberGen.m_Translation[2]); m_FiberGen.m_Scale[0] = ReadVal(v1,"scale.x", m_FiberGen.m_Scale[0]); m_FiberGen.m_Scale[1] = ReadVal(v1,"scale.y", m_FiberGen.m_Scale[1]); m_FiberGen.m_Scale[2] = ReadVal(v1,"scale.z", m_FiberGen.m_Scale[2]); } else if ( v1.first == "image" ) { m_Misc.m_SignalModelString = ReadVal(v1,"signalmodelstring", m_Misc.m_SignalModelString); m_Misc.m_ArtifactModelString = ReadVal(v1,"artifactmodelstring", m_Misc.m_ArtifactModelString); m_Misc.m_OutputPath = ReadVal(v1,"outpath", m_Misc.m_OutputPath); m_Misc.m_CheckOutputVolumeFractionsBox = ReadVal(v1,"outputvolumefractions", m_Misc.m_CheckOutputVolumeFractionsBox); m_Misc.m_CheckAdvancedSignalOptionsBox = ReadVal(v1,"showadvanced", m_Misc.m_CheckAdvancedSignalOptionsBox); m_Misc.m_DoAddDistortions = ReadVal(v1,"artifacts.doAddDistortions", m_Misc.m_DoAddDistortions); m_Misc.m_DoAddNoise = ReadVal(v1,"artifacts.addnoise", m_Misc.m_DoAddNoise); m_Misc.m_DoAddGhosts = ReadVal(v1,"artifacts.addghosts", m_Misc.m_DoAddGhosts); m_Misc.m_DoAddAliasing = ReadVal(v1,"artifacts.addaliasing", m_Misc.m_DoAddAliasing); m_Misc.m_DoAddSpikes = ReadVal(v1,"artifacts.addspikes", m_Misc.m_DoAddSpikes); m_Misc.m_DoAddEddyCurrents = ReadVal(v1,"artifacts.addeddycurrents", m_Misc.m_DoAddEddyCurrents); m_SignalGen.m_ImageRegion.SetSize(0, ReadVal(v1,"basic.size.x",m_SignalGen.m_ImageRegion.GetSize(0))); m_SignalGen.m_ImageRegion.SetSize(1, ReadVal(v1,"basic.size.y",m_SignalGen.m_ImageRegion.GetSize(1))); m_SignalGen.m_ImageRegion.SetSize(2, ReadVal(v1,"basic.size.z",m_SignalGen.m_ImageRegion.GetSize(2))); m_SignalGen.m_ImageSpacing[0] = ReadVal(v1,"basic.spacing.x",m_SignalGen.m_ImageSpacing[0]); m_SignalGen.m_ImageSpacing[1] = ReadVal(v1,"basic.spacing.y",m_SignalGen.m_ImageSpacing[1]); m_SignalGen.m_ImageSpacing[2] = ReadVal(v1,"basic.spacing.z",m_SignalGen.m_ImageSpacing[2]); m_SignalGen.m_ImageOrigin[0] = ReadVal(v1,"basic.origin.x",m_SignalGen.m_ImageOrigin[0]); m_SignalGen.m_ImageOrigin[1] = ReadVal(v1,"basic.origin.y",m_SignalGen.m_ImageOrigin[1]); m_SignalGen.m_ImageOrigin[2] = ReadVal(v1,"basic.origin.z",m_SignalGen.m_ImageOrigin[2]); int i = 0; int j = 0; for(auto v : v1.second.get_child("basic.direction")) { m_SignalGen.m_ImageDirection[i][j] = boost::lexical_cast(v.second.data()); ++j; if (j==3) { j = 0; ++i; } } m_SignalGen.m_AcquisitionType = (SignalGenerationParameters::AcquisitionType)ReadVal(v1,"acquisitiontype", m_SignalGen.m_AcquisitionType); m_SignalGen.m_CoilSensitivityProfile = (SignalGenerationParameters::CoilSensitivityProfile)ReadVal(v1,"coilsensitivityprofile", m_SignalGen.m_CoilSensitivityProfile); m_SignalGen.m_NumberOfCoils = ReadVal(v1,"numberofcoils", m_SignalGen.m_NumberOfCoils); m_SignalGen.m_ReversePhase = ReadVal(v1,"reversephase", m_SignalGen.m_ReversePhase); m_SignalGen.m_PartialFourier = ReadVal(v1,"partialfourier", m_SignalGen.m_PartialFourier); m_SignalGen.m_NoiseVariance = ReadVal(v1,"noisevariance", m_SignalGen.m_NoiseVariance); m_SignalGen.m_tRep = ReadVal(v1,"trep", m_SignalGen.m_tRep); m_SignalGen.m_SignalScale = ReadVal(v1,"signalScale", m_SignalGen.m_SignalScale); m_SignalGen.m_tEcho = ReadVal(v1,"tEcho", m_SignalGen.m_tEcho); m_SignalGen.m_tLine = ReadVal(v1,"tLine", m_SignalGen.m_tLine); m_SignalGen.m_tInhom = ReadVal(v1,"tInhom", m_SignalGen.m_tInhom); m_SignalGen.m_SimulateKspaceAcquisition = ReadVal(v1,"simulatekspace", m_SignalGen.m_SimulateKspaceAcquisition); m_SignalGen.m_AxonRadius = ReadVal(v1,"axonRadius", m_SignalGen.m_AxonRadius); m_SignalGen.m_Spikes = ReadVal(v1,"artifacts.spikesnum", m_SignalGen.m_Spikes); m_SignalGen.m_SpikeAmplitude = ReadVal(v1,"artifacts.spikesscale", m_SignalGen.m_SpikeAmplitude); m_SignalGen.m_KspaceLineOffset = ReadVal(v1,"artifacts.kspaceLineOffset", m_SignalGen.m_KspaceLineOffset); m_SignalGen.m_EddyStrength = ReadVal(v1,"artifacts.eddyStrength", m_SignalGen.m_EddyStrength); m_SignalGen.m_Tau = ReadVal(v1,"artifacts.eddyTau", m_SignalGen.m_Tau); m_SignalGen.m_CroppingFactor = ReadVal(v1,"artifacts.aliasingfactor", m_SignalGen.m_CroppingFactor); m_SignalGen.m_Drift = ReadVal(v1,"artifacts.drift", m_SignalGen.m_Drift); m_SignalGen.m_DoAddGibbsRinging = ReadVal(v1,"artifacts.addringing", m_SignalGen.m_DoAddGibbsRinging); + m_SignalGen.m_ZeroRinging = ReadVal(v1,"artifacts.zeroringing", m_SignalGen.m_ZeroRinging); m_SignalGen.m_DoSimulateRelaxation = ReadVal(v1,"doSimulateRelaxation", m_SignalGen.m_DoSimulateRelaxation); m_SignalGen.m_DoDisablePartialVolume = ReadVal(v1,"doDisablePartialVolume", m_SignalGen.m_DoDisablePartialVolume); m_SignalGen.m_DoAddMotion = ReadVal(v1,"artifacts.doAddMotion", m_SignalGen.m_DoAddMotion); m_SignalGen.m_DoRandomizeMotion = ReadVal(v1,"artifacts.randomMotion", m_SignalGen.m_DoRandomizeMotion); m_SignalGen.m_DoAddDrift = ReadVal(v1,"artifacts.doAddDrift", m_SignalGen.m_DoAddDrift); m_SignalGen.m_Translation[0] = ReadVal(v1,"artifacts.translation0", m_SignalGen.m_Translation[0]); m_SignalGen.m_Translation[1] = ReadVal(v1,"artifacts.translation1", m_SignalGen.m_Translation[1]); m_SignalGen.m_Translation[2] = ReadVal(v1,"artifacts.translation2", m_SignalGen.m_Translation[2]); m_SignalGen.m_Rotation[0] = ReadVal(v1,"artifacts.rotation0", m_SignalGen.m_Rotation[0]); m_SignalGen.m_Rotation[1] = ReadVal(v1,"artifacts.rotation1", m_SignalGen.m_Rotation[1]); m_SignalGen.m_Rotation[2] = ReadVal(v1,"artifacts.rotation2", m_SignalGen.m_Rotation[2]); if (itksys::SystemTools::FileExists(filename+".bvals") && itksys::SystemTools::FileExists(filename+".bvecs")) { m_Misc.m_BvalsFile = filename+".bvals"; m_Misc.m_BvecsFile = filename+".bvecs"; m_SignalGen.SetGradienDirections( mitk::gradients::ReadBvalsBvecs(m_Misc.m_BvalsFile, m_Misc.m_BvecsFile, m_SignalGen.m_Bvalue) ); } else { m_SignalGen.m_Bvalue = ReadVal(v1,"bvalue", m_SignalGen.m_Bvalue); SignalGenerationParameters::GradientListType gradients; try { BOOST_FOREACH( boost::property_tree::ptree::value_type const& v2, v1.second.get_child("gradients") ) { SignalGenerationParameters::GradientType g; g[0] = ReadVal(v2,"x",0); g[1] = ReadVal(v2,"y",0); g[2] = ReadVal(v2,"z",0); gradients.push_back(g); } } catch(...) { MITK_INFO << "WARNING: Fiberfox parameters without any gradient directions loaded."; } m_SignalGen.SetGradienDirections(gradients); } m_Misc.m_MotionVolumesBox = ReadVal(v1,"artifacts.motionvolumes", m_Misc.m_MotionVolumesBox); m_SignalGen.m_MotionVolumes.clear(); if ( m_Misc.m_MotionVolumesBox == "random" ) { m_SignalGen.m_MotionVolumes.push_back(0); for ( size_t i=1; i < m_SignalGen.GetNumVolumes(); ++i ) { m_SignalGen.m_MotionVolumes.push_back( bool( randgen->GetIntegerVariate()%2 ) ); } MITK_DEBUG << "mitkFiberfoxParameters.cpp: Case m_Misc.m_MotionVolumesBox == \"random\"."; } else if ( ! m_Misc.m_MotionVolumesBox.empty() ) { std::stringstream stream( m_Misc.m_MotionVolumesBox ); std::vector numbers; int nummer = std::numeric_limits::max(); while( stream >> nummer ) { if( nummer < std::numeric_limits::max() ) { numbers.push_back( nummer ); } } // If a list of negative numbers is given: if( *(std::min_element( numbers.begin(), numbers.end() )) < 0 && *(std::max_element( numbers.begin(), numbers.end() )) <= 0 ) // cave: -0 == +0 { for ( size_t i=0; i(m_SignalGen.GetNumVolumes()) && -number >= 0 ) m_SignalGen.m_MotionVolumes.at(-number) = false; } MITK_DEBUG << "mitkFiberfoxParameters.cpp: Case list of negative numbers."; } // If a list of positive numbers is given: else if( *(std::min_element( numbers.begin(), numbers.end() )) >= 0 && *(std::max_element( numbers.begin(), numbers.end() )) >= 0 ) { for ( size_t i=0; i(m_SignalGen.GetNumVolumes()) && number >= 0) m_SignalGen.m_MotionVolumes.at(number) = true; } MITK_DEBUG << "mitkFiberfoxParameters.cpp: Case list of positive numbers."; } else { MITK_WARN << "mitkFiberfoxParameters.cpp: Inconsistent list of numbers in m_MotionVolumesBox."; break; } } else { MITK_WARN << "mitkFiberfoxParameters.cpp: Cannot make sense of string in m_MotionVolumesBox."; break; } try { if (ReadVal(v1,"artifacts.noisetype","")=="rice") { m_NoiseModel = std::make_shared< mitk::RicianNoiseModel<> >(); m_NoiseModel->SetNoiseVariance(ReadVal(v1,"artifacts.noisevariance",m_NoiseModel->GetNoiseVariance())); } } catch(...) { MITK_DEBUG << "mitkFiberfoxParameters.cpp: caught some error while trying m_NoiseModel->SetNoiseVariance()"; // throw; } try { if (ReadVal(v1,"artifacts.noisetype","")=="chisquare") { m_NoiseModel = std::make_shared< mitk::ChiSquareNoiseModel<> >(); m_NoiseModel->SetNoiseVariance(ReadVal(v1,"artifacts.noisevariance",m_NoiseModel->GetNoiseVariance())); } } catch(...) { MITK_DEBUG << "mitkFiberfoxParameters.cpp: caught some error while trying m_NoiseModel->SetNoiseVariance()"; // throw; } BOOST_FOREACH( boost::property_tree::ptree::value_type const& v2, v1.second.get_child("compartments") ) { mitk::DiffusionSignalModel<>* signalModel = nullptr; std::string model = ReadVal(v2,"model","",true); if (model=="stick") { mitk::StickModel<>* model = new mitk::StickModel<>(); model->SetDiffusivity(ReadVal(v2,"d",model->GetDiffusivity())); model->SetT2(ReadVal(v2,"t2",model->GetT2())); model->SetT1(ReadVal(v2,"t1",model->GetT1())); model->SetBvalue(m_SignalGen.m_Bvalue); model->m_CompartmentId = ReadVal(v2,"ID",0,true); if (ReadVal(v2,"type","",true)=="fiber") m_FiberModelList.push_back(model); else if (ReadVal(v2,"type","",true)=="non-fiber") m_NonFiberModelList.push_back(model); signalModel = model; } else if (model=="tensor") { mitk::TensorModel<>* model = new mitk::TensorModel<>(); model->SetDiffusivity1(ReadVal(v2,"d1",model->GetDiffusivity1())); model->SetDiffusivity2(ReadVal(v2,"d2",model->GetDiffusivity2())); model->SetDiffusivity3(ReadVal(v2,"d3",model->GetDiffusivity3())); model->SetT2(ReadVal(v2,"t2",model->GetT2())); model->SetT1(ReadVal(v2,"t1",model->GetT1())); model->SetBvalue(m_SignalGen.m_Bvalue); model->m_CompartmentId = ReadVal(v2,"ID",0,true); if (ReadVal(v2,"type","",true)=="fiber") m_FiberModelList.push_back(model); else if (ReadVal(v2,"type","",true)=="non-fiber") m_NonFiberModelList.push_back(model); signalModel = model; } else if (model=="ball") { mitk::BallModel<>* model = new mitk::BallModel<>(); model->SetDiffusivity(ReadVal(v2,"d",model->GetDiffusivity())); model->SetT2(ReadVal(v2,"t2",model->GetT2())); model->SetT1(ReadVal(v2,"t1",model->GetT1())); model->SetBvalue(m_SignalGen.m_Bvalue); model->m_CompartmentId = ReadVal(v2,"ID",0,true); if (ReadVal(v2,"type","",true)=="fiber") m_FiberModelList.push_back(model); else if (ReadVal(v2,"type","",true)=="non-fiber") m_NonFiberModelList.push_back(model); signalModel = model; } else if (model=="astrosticks") { mitk::AstroStickModel<>* model = new AstroStickModel<>(); model->SetDiffusivity(ReadVal(v2,"d",model->GetDiffusivity())); model->SetT2(ReadVal(v2,"t2",model->GetT2())); model->SetT1(ReadVal(v2,"t1",model->GetT1())); model->SetBvalue(m_SignalGen.m_Bvalue); model->SetRandomizeSticks(ReadVal(v2,"randomize",model->GetRandomizeSticks())); model->m_CompartmentId = ReadVal(v2,"ID",0,true); if (ReadVal(v2,"type","",true)=="fiber") m_FiberModelList.push_back(model); else if (ReadVal(v2,"type","",true)=="non-fiber") m_NonFiberModelList.push_back(model); signalModel = model; } else if (model=="dot") { mitk::DotModel<>* model = new mitk::DotModel<>(); model->SetT2(ReadVal(v2,"t2",model->GetT2())); model->SetT1(ReadVal(v2,"t1",model->GetT1())); model->m_CompartmentId = ReadVal(v2,"ID",0,true); if (ReadVal(v2,"type","",true)=="fiber") m_FiberModelList.push_back(model); else if (ReadVal(v2,"type","",true)=="non-fiber") m_NonFiberModelList.push_back(model); signalModel = model; } else if (model=="prototype") { mitk::RawShModel<>* model = new mitk::RawShModel<>(); model->SetMaxNumKernels(ReadVal(v2,"maxNumSamples",model->GetMaxNumKernels())); model->SetFaRange(ReadVal(v2,"minFA",model->GetFaRange().first), ReadVal(v2,"maxFA",model->GetFaRange().second)); model->SetAdcRange(ReadVal(v2,"minADC",model->GetAdcRange().first), ReadVal(v2,"maxADC",model->GetAdcRange().second)); model->m_CompartmentId = ReadVal(v2,"ID",0,true); unsigned int numCoeffs = ReadVal(v2,"numCoeffs",0,true); unsigned int numSamples = ReadVal(v2,"numSamples",0,true); for (unsigned int j=0; j coeffs(numCoeffs); for (unsigned int k=0; k(v2,"kernels."+boost::lexical_cast(j)+".coeffs."+boost::lexical_cast(k),0,true); } model->SetShCoefficients( coeffs, ReadVal(v2,"kernels."+boost::lexical_cast(j)+".B0",0,true) ); } if (ReadVal(v2,"type","",true)=="fiber") { m_FiberModelList.push_back(model); } else if (ReadVal(v2,"type","",true)=="non-fiber") { m_NonFiberModelList.push_back(model); } // else ? signalModel = model; } if (signalModel!=nullptr) { try { itk::ImageFileReader::Pointer reader = itk::ImageFileReader::New(); if ( itksys::SystemTools::FileExists(filename+"_VOLUME"+ReadVal(v2,"ID","")+".nii.gz") ) reader->SetFileName(filename+"_VOLUME"+ReadVal(v2,"ID","")+".nii.gz"); else if ( itksys::SystemTools::FileExists(filename+"_VOLUME"+ReadVal(v2,"ID","")+".nii") ) reader->SetFileName(filename+"_VOLUME"+ReadVal(v2,"ID","")+".nii"); else reader->SetFileName(filename+"_VOLUME"+ReadVal(v2,"ID","")+".nrrd"); reader->Update(); signalModel->SetVolumeFractionImage(reader->GetOutput()); MITK_INFO << "Volume fraction image loaded for compartment " << signalModel->m_CompartmentId; } catch(...) { MITK_INFO << "No volume fraction image found for compartment " << signalModel->m_CompartmentId; } } } } else { } } UpdateSignalModels(); try { itk::ImageFileReader::Pointer reader = itk::ImageFileReader::New(); reader->SetFileName(filename+"_FMAP.nrrd"); if ( itksys::SystemTools::FileExists(filename+"_FMAP.nii.gz") ) reader->SetFileName(filename+"_FMAP.nii.gz"); else if ( itksys::SystemTools::FileExists(filename+"_FMAP.nii") ) reader->SetFileName(filename+"_FMAP.nii"); else reader->SetFileName(filename+"_FMAP.nrrd"); reader->Update(); m_SignalGen.m_FrequencyMap = reader->GetOutput(); MITK_INFO << "Frequency map loaded."; } catch(...) { MITK_INFO << "No frequency map found."; } try { itk::ImageFileReader::Pointer reader = itk::ImageFileReader::New(); if ( itksys::SystemTools::FileExists(filename+"_MASK.nii.gz") ) reader->SetFileName(filename+"_MASK.nii.gz"); else if ( itksys::SystemTools::FileExists(filename+"_MASK.nii") ) reader->SetFileName(filename+"_MASK.nii"); else reader->SetFileName(filename+"_MASK.nrrd"); reader->Update(); m_SignalGen.m_MaskImage = reader->GetOutput(); m_SignalGen.m_ImageRegion = m_SignalGen.m_MaskImage->GetLargestPossibleRegion(); m_SignalGen.m_ImageSpacing = m_SignalGen.m_MaskImage->GetSpacing(); m_SignalGen.m_ImageOrigin = m_SignalGen.m_MaskImage->GetOrigin(); m_SignalGen.m_ImageDirection = m_SignalGen.m_MaskImage->GetDirection(); MITK_INFO << "Mask image loaded."; } catch(...) { MITK_INFO << "No mask image found."; } setlocale(LC_ALL, currLocale.c_str()); } void mitk::FiberfoxParameters::PrintSelf() { MITK_INFO << "Not implemented :("; } diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h index 72304a9356..c219688526 100644 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h @@ -1,329 +1,331 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _MITK_FiberfoxParameters_H #define _MITK_FiberfoxParameters_H #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #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 }; SignalGenerationParameters() : m_AcquisitionType(SignalGenerationParameters::SingleShotEpi) , m_SignalScale(100) , m_tEcho(100) , m_tRep(4000) , 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_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 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/Plugins/org.mitk.gui.qt.diffusionimaging.fiberfox/src/internal/QmitkFiberfoxView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberfox/src/internal/QmitkFiberfoxView.cpp index 592fe77b59..60ff07e39c 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberfox/src/internal/QmitkFiberfoxView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberfox/src/internal/QmitkFiberfoxView.cpp @@ -1,2127 +1,2139 @@ /*=================================================================== 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. ===================================================================*/ // Blueberry #include #include // Qmitk #include "QmitkFiberfoxView.h" // MITK #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #define RAPIDXML_NO_EXCEPTIONS #include #include #include #include #include "usModuleRegistry.h" #include #include #include #include #include #include #include #include #include #include "mitkNodePredicateDataType.h" #include #include #include #include #define _USE_MATH_DEFINES #include QmitkFiberfoxWorker::QmitkFiberfoxWorker(QmitkFiberfoxView* view) : m_View(view) { } void QmitkFiberfoxWorker::run() { try{ m_View->m_TractsToDwiFilter->Update(); } catch( ... ) { } m_View->m_Thread.quit(); } const std::string QmitkFiberfoxView::VIEW_ID = "org.mitk.views.fiberfoxview"; QmitkFiberfoxView::QmitkFiberfoxView() : QmitkAbstractView() , m_Controls( 0 ) , m_SelectedImageNode( nullptr ) , m_Worker(this) , m_ThreadIsRunning(false) { m_Worker.moveToThread(&m_Thread); connect(&m_Thread, SIGNAL(started()), this, SLOT(BeforeThread())); connect(&m_Thread, SIGNAL(started()), &m_Worker, SLOT(run())); connect(&m_Thread, SIGNAL(finished()), this, SLOT(AfterThread())); // connect(&m_Thread, SIGNAL(terminated()), this, SLOT(AfterThread())); m_SimulationTimer = new QTimer(this); } void QmitkFiberfoxView::KillThread() { MITK_INFO << "Aborting DWI simulation."; m_TractsToDwiFilter->SetAbortGenerateData(true); m_Controls->m_AbortSimulationButton->setEnabled(false); m_Controls->m_AbortSimulationButton->setText("Aborting simulation ..."); } void QmitkFiberfoxView::BeforeThread() { m_SimulationTime = QTime::currentTime(); m_SimulationTimer->start(100); m_Controls->m_AbortSimulationButton->setVisible(true); m_Controls->m_GenerateImageButton->setVisible(false); m_Controls->m_SimulationStatusText->setVisible(true); m_ThreadIsRunning = true; } void QmitkFiberfoxView::AfterThread() { UpdateSimulationStatus(); m_SimulationTimer->stop(); m_Controls->m_AbortSimulationButton->setVisible(false); m_Controls->m_AbortSimulationButton->setEnabled(true); m_Controls->m_AbortSimulationButton->setText("Abort simulation"); m_Controls->m_GenerateImageButton->setVisible(true); m_ThreadIsRunning = false; QString statusText; FiberfoxParameters parameters; mitk::Image::Pointer mitkImage = mitk::Image::New(); statusText = QString(m_TractsToDwiFilter->GetStatusText().c_str()); if (m_TractsToDwiFilter->GetAbortGenerateData()) { MITK_INFO << "Simulation aborted."; return; } parameters = m_TractsToDwiFilter->GetParameters(); mitkImage = mitk::GrabItkImageMemory( m_TractsToDwiFilter->GetOutput() ); mitk::DiffusionPropertyHelper::SetGradientContainer(mitkImage, parameters.m_SignalGen.GetItkGradientContainer()); mitk::DiffusionPropertyHelper::SetReferenceBValue(mitkImage, parameters.m_SignalGen.GetBvalue()); mitk::DiffusionPropertyHelper::InitializeImage( mitkImage ); parameters.m_Misc.m_ResultNode->SetData( mitkImage ); GetDataStorage()->Add(parameters.m_Misc.m_ResultNode, parameters.m_Misc.m_ParentNode); parameters.m_Misc.m_ResultNode->SetProperty( "levelwindow", mitk::LevelWindowProperty::New(m_TractsToDwiFilter->GetLevelWindow()) ); if (m_Controls->m_VolumeFractionsBox->isChecked()) { if (m_TractsToDwiFilter->GetPhaseImage().IsNotNull()) { mitk::Image::Pointer phaseImage = mitk::Image::New(); itk::TractsToDWIImageFilter< short >::DoubleDwiType::Pointer itkPhase = m_TractsToDwiFilter->GetPhaseImage(); phaseImage = mitk::GrabItkImageMemory( itkPhase.GetPointer() ); mitk::DataNode::Pointer phaseNode = mitk::DataNode::New(); phaseNode->SetData( phaseImage ); phaseNode->SetName("Phase Image"); GetDataStorage()->Add(phaseNode, parameters.m_Misc.m_ResultNode); } if (m_TractsToDwiFilter->GetKspaceImage().IsNotNull()) { mitk::Image::Pointer image = mitk::Image::New(); itk::TractsToDWIImageFilter< short >::DoubleDwiType::Pointer itkImage = m_TractsToDwiFilter->GetKspaceImage(); image = mitk::GrabItkImageMemory( itkImage.GetPointer() ); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); node->SetName("k-Space"); GetDataStorage()->Add(node, parameters.m_Misc.m_ResultNode); } { mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(m_TractsToDwiFilter->GetCoilPointset()); node->SetName("Coil Positions"); node->SetProperty("pointsize", mitk::FloatProperty::New(parameters.m_SignalGen.m_ImageSpacing[0]/4)); node->SetProperty("color", mitk::ColorProperty::New(0, 1, 0)); GetDataStorage()->Add(node, parameters.m_Misc.m_ResultNode); } int c = 1; std::vector< itk::TractsToDWIImageFilter< short >::DoubleDwiType::Pointer > output_real = m_TractsToDwiFilter->GetOutputImagesReal(); for (auto real : output_real) { mitk::Image::Pointer image = mitk::Image::New(); image->InitializeByItk(real.GetPointer()); image->SetVolume(real->GetBufferPointer()); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); node->SetName("Coil-"+QString::number(c).toStdString()+"-real"); GetDataStorage()->Add(node, parameters.m_Misc.m_ResultNode); ++c; } c = 1; std::vector< itk::TractsToDWIImageFilter< short >::DoubleDwiType::Pointer > output_imag = m_TractsToDwiFilter->GetOutputImagesImag(); for (auto imag : output_imag) { mitk::Image::Pointer image = mitk::Image::New(); image->InitializeByItk(imag.GetPointer()); image->SetVolume(imag->GetBufferPointer()); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); node->SetName("Coil-"+QString::number(c).toStdString()+"-imag"); GetDataStorage()->Add(node, parameters.m_Misc.m_ResultNode); ++c; } std::vector< itk::TractsToDWIImageFilter< short >::ItkDoubleImgType::Pointer > volumeFractions = m_TractsToDwiFilter->GetVolumeFractions(); for (unsigned int k=0; kInitializeByItk(volumeFractions.at(k).GetPointer()); image->SetVolume(volumeFractions.at(k)->GetBufferPointer()); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); node->SetName("CompartmentVolume-"+QString::number(k).toStdString()); GetDataStorage()->Add(node, parameters.m_Misc.m_ResultNode); } } m_TractsToDwiFilter = nullptr; if (parameters.m_Misc.m_AfterSimulationMessage.size()>0) QMessageBox::information( nullptr, "Warning", parameters.m_Misc.m_AfterSimulationMessage.c_str()); mitk::BaseData::Pointer basedata = parameters.m_Misc.m_ResultNode->GetData(); if (basedata.IsNotNull()) { mitk::RenderingManager::GetInstance()->InitializeViews( basedata->GetTimeGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true ); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } if (!parameters.m_Misc.m_OutputPath.empty()) { try{ QString outputFileName(parameters.m_Misc.m_OutputPath.c_str()); outputFileName += parameters.m_Misc.m_ResultNode->GetName().c_str(); outputFileName.replace(QString("."), QString("_")); SaveParameters(outputFileName+".ffp"); outputFileName += ".dwi"; QString status("Saving output image to "); status += outputFileName; m_Controls->m_SimulationStatusText->append(status); mitk::IOUtil::Save(mitkImage, outputFileName.toStdString()); m_Controls->m_SimulationStatusText->append("File saved successfully."); } catch (itk::ExceptionObject &e) { QString status("Exception during DWI writing: "); status += e.GetDescription(); m_Controls->m_SimulationStatusText->append(status); } catch (...) { m_Controls->m_SimulationStatusText->append("Unknown exception during DWI writing!"); } } parameters.m_SignalGen.m_FrequencyMap = nullptr; } void QmitkFiberfoxView::UpdateSimulationStatus() { QString statusText = QString(m_TractsToDwiFilter->GetStatusText().c_str()); if (QString::compare(m_SimulationStatusText,statusText)!=0) { m_Controls->m_SimulationStatusText->clear(); m_Controls->m_SimulationStatusText->setText(statusText); QScrollBar *vScrollBar = m_Controls->m_SimulationStatusText->verticalScrollBar(); vScrollBar->triggerAction(QScrollBar::SliderToMaximum); } } // Destructor QmitkFiberfoxView::~QmitkFiberfoxView() { delete m_SimulationTimer; } void QmitkFiberfoxView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkFiberfoxViewControls; m_Controls->setupUi( parent ); m_Controls->m_StickWidget1->setVisible(true); m_Controls->m_StickWidget2->setVisible(false); m_Controls->m_ZeppelinWidget1->setVisible(false); m_Controls->m_ZeppelinWidget2->setVisible(false); m_Controls->m_TensorWidget1->setVisible(false); m_Controls->m_TensorWidget2->setVisible(false); m_Controls->m_BallWidget1->setVisible(true); m_Controls->m_BallWidget2->setVisible(false); m_Controls->m_BallWidget2->SetT1(4500); m_Controls->m_AstrosticksWidget1->setVisible(false); m_Controls->m_AstrosticksWidget2->setVisible(false); m_Controls->m_AstrosticksWidget2->SetT1(4500); m_Controls->m_DotWidget1->setVisible(false); m_Controls->m_DotWidget2->setVisible(false); m_Controls->m_DotWidget2->SetT1(4500); m_Controls->m_PrototypeWidget1->setVisible(false); m_Controls->m_PrototypeWidget2->setVisible(false); m_Controls->m_PrototypeWidget3->setVisible(false); m_Controls->m_PrototypeWidget4->setVisible(false); m_Controls->m_PrototypeWidget3->SetMinFa(0.0); m_Controls->m_PrototypeWidget3->SetMaxFa(0.15); m_Controls->m_PrototypeWidget4->SetMinFa(0.0); m_Controls->m_PrototypeWidget4->SetMaxFa(0.15); m_Controls->m_PrototypeWidget3->SetMinAdc(0.0); m_Controls->m_PrototypeWidget3->SetMaxAdc(0.001); m_Controls->m_PrototypeWidget4->SetMinAdc(0.003); m_Controls->m_PrototypeWidget4->SetMaxAdc(0.004); m_Controls->m_Comp2FractionFrame->setVisible(false); m_Controls->m_Comp4FractionFrame->setVisible(false); m_Controls->m_DiffusionPropsMessage->setVisible(false); m_Controls->m_GeometryMessage->setVisible(false); m_Controls->m_AdvancedSignalOptionsFrame->setVisible(false); m_Controls->m_NoiseFrame->setVisible(false); + m_Controls->m_ZeroRinging->setVisible(false); m_Controls->m_GhostFrame->setVisible(false); m_Controls->m_DistortionsFrame->setVisible(false); m_Controls->m_EddyFrame->setVisible(false); m_Controls->m_SpikeFrame->setVisible(false); m_Controls->m_AliasingFrame->setVisible(false); m_Controls->m_MotionArtifactFrame->setVisible(false); m_Controls->m_DriftFrame->setVisible(false); m_ParameterFile = QDir::currentPath()+"/param.ffp"; m_Controls->m_AbortSimulationButton->setVisible(false); m_Controls->m_SimulationStatusText->setVisible(false); m_Controls->m_FrequencyMapBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_Comp1VolumeFraction->SetDataStorage(this->GetDataStorage()); m_Controls->m_Comp2VolumeFraction->SetDataStorage(this->GetDataStorage()); m_Controls->m_Comp3VolumeFraction->SetDataStorage(this->GetDataStorage()); m_Controls->m_Comp4VolumeFraction->SetDataStorage(this->GetDataStorage()); m_Controls->m_MaskComboBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_TemplateComboBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_FiberBundleComboBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isFiberBundle = mitk::TNodePredicateDataType::New(); mitk::TNodePredicateDataType::Pointer isMitkImage = mitk::TNodePredicateDataType::New(); mitk::NodePredicateIsDWI::Pointer isDwi = mitk::NodePredicateIsDWI::New( ); mitk::NodePredicateDataType::Pointer isDti = mitk::NodePredicateDataType::New("TensorImage"); mitk::NodePredicateDataType::Pointer isOdf = mitk::NodePredicateDataType::New("Odfmage"); mitk::NodePredicateOr::Pointer isDiffusionImage = mitk::NodePredicateOr::New(isDwi, isDti); isDiffusionImage = mitk::NodePredicateOr::New(isDiffusionImage, isOdf); mitk::NodePredicateNot::Pointer noDiffusionImage = mitk::NodePredicateNot::New(isDiffusionImage); mitk::NodePredicateAnd::Pointer isNonDiffMitkImage = mitk::NodePredicateAnd::New(isMitkImage, noDiffusionImage); mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); mitk::NodePredicateAnd::Pointer isBinaryMitkImage = mitk::NodePredicateAnd::New( isNonDiffMitkImage, isBinaryPredicate ); m_Controls->m_FrequencyMapBox->SetPredicate(isNonDiffMitkImage); m_Controls->m_Comp1VolumeFraction->SetPredicate(isNonDiffMitkImage); m_Controls->m_Comp1VolumeFraction->SetZeroEntryText("--"); m_Controls->m_Comp2VolumeFraction->SetPredicate(isNonDiffMitkImage); m_Controls->m_Comp2VolumeFraction->SetZeroEntryText("--"); m_Controls->m_Comp3VolumeFraction->SetPredicate(isNonDiffMitkImage); m_Controls->m_Comp3VolumeFraction->SetZeroEntryText("--"); m_Controls->m_Comp4VolumeFraction->SetPredicate(isNonDiffMitkImage); m_Controls->m_Comp4VolumeFraction->SetZeroEntryText("--"); m_Controls->m_MaskComboBox->SetPredicate(isBinaryMitkImage); m_Controls->m_MaskComboBox->SetZeroEntryText("--"); m_Controls->m_TemplateComboBox->SetPredicate(isMitkImage); m_Controls->m_TemplateComboBox->SetZeroEntryText("--"); m_Controls->m_FiberBundleComboBox->SetPredicate(isFiberBundle); m_Controls->m_FiberBundleComboBox->SetZeroEntryText("--"); QFont font; font.setFamily("Courier"); font.setStyleHint(QFont::Monospace); font.setFixedPitch(true); font.setPointSize(7); m_Controls->m_SimulationStatusText->setFont(font); connect( m_SimulationTimer, SIGNAL(timeout()), this, SLOT(UpdateSimulationStatus()) ); connect((QObject*) m_Controls->m_AbortSimulationButton, SIGNAL(clicked()), (QObject*) this, SLOT(KillThread())); connect((QObject*) m_Controls->m_GenerateImageButton, SIGNAL(clicked()), (QObject*) this, SLOT(GenerateImage())); connect((QObject*) m_Controls->m_AddNoise, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddNoise(int))); connect((QObject*) m_Controls->m_AddGhosts, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddGhosts(int))); connect((QObject*) m_Controls->m_AddDistortions, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddDistortions(int))); connect((QObject*) m_Controls->m_AddEddy, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddEddy(int))); connect((QObject*) m_Controls->m_AddSpikes, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddSpikes(int))); connect((QObject*) m_Controls->m_AddAliasing, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddAliasing(int))); connect((QObject*) m_Controls->m_AddMotion, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddMotion(int))); connect((QObject*) m_Controls->m_AddDrift, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddDrift(int))); + connect((QObject*) m_Controls->m_AddGibbsRinging, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddRinging(int))); connect((QObject*) m_Controls->m_Compartment1Box, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(Comp1ModelFrameVisibility(int))); connect((QObject*) m_Controls->m_Compartment2Box, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(Comp2ModelFrameVisibility(int))); connect((QObject*) m_Controls->m_Compartment3Box, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(Comp3ModelFrameVisibility(int))); connect((QObject*) m_Controls->m_Compartment4Box, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(Comp4ModelFrameVisibility(int))); connect((QObject*) m_Controls->m_AdvancedOptionsBox_2, SIGNAL( stateChanged(int)), (QObject*) this, SLOT(ShowAdvancedOptions(int))); connect((QObject*) m_Controls->m_UseBvalsBvecsBox, SIGNAL( stateChanged(int)), (QObject*) this, SLOT(OnBvalsBvecsCheck(int))); connect((QObject*) m_Controls->m_SaveParametersButton, SIGNAL(clicked()), (QObject*) this, SLOT(SaveParameters())); connect((QObject*) m_Controls->m_LoadParametersButton, SIGNAL(clicked()), (QObject*) this, SLOT(LoadParameters())); connect((QObject*) m_Controls->m_OutputPathButton, SIGNAL(clicked()), (QObject*) this, SLOT(SetOutputPath())); connect((QObject*) m_Controls->m_LoadBvalsButton, SIGNAL(clicked()), (QObject*) this, SLOT(SetBvalsEdit())); connect((QObject*) m_Controls->m_LoadBvecsButton, SIGNAL(clicked()), (QObject*) this, SLOT(SetBvecsEdit())); connect((QObject*) m_Controls->m_MaskComboBox, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(OnMaskSelected(int))); connect((QObject*) m_Controls->m_TemplateComboBox, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(OnTemplateSelected(int))); connect((QObject*) m_Controls->m_FiberBundleComboBox, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(OnFibSelected(int))); } UpdateGui(); } void QmitkFiberfoxView::OnMaskSelected(int ) { UpdateGui(); } void QmitkFiberfoxView::OnTemplateSelected(int ) { UpdateGui(); } void QmitkFiberfoxView::OnFibSelected(int ) { UpdateGui(); } void QmitkFiberfoxView::OnBvalsBvecsCheck(int ) { UpdateGui(); } void QmitkFiberfoxView::UpdateParametersFromGui() { m_Parameters.ClearSignalParameters(); m_Parameters.m_Misc.m_CheckAdvancedSignalOptionsBox = m_Controls->m_AdvancedOptionsBox_2->isChecked(); m_Parameters.m_Misc.m_CheckOutputVolumeFractionsBox = m_Controls->m_VolumeFractionsBox->isChecked(); std::string outputPath = m_Controls->m_SavePathEdit->text().toStdString(); if (outputPath.compare("-")!=0) { m_Parameters.m_Misc.m_OutputPath = outputPath; m_Parameters.m_Misc.m_OutputPath += "/"; } if (m_Controls->m_MaskComboBox->GetSelectedNode().IsNotNull()) { mitk::Image::Pointer mitkMaskImage = dynamic_cast(m_Controls->m_MaskComboBox->GetSelectedNode()->GetData()); mitk::CastToItkImage(mitkMaskImage, m_Parameters.m_SignalGen.m_MaskImage); itk::ImageDuplicator::Pointer duplicator = itk::ImageDuplicator::New(); duplicator->SetInputImage(m_Parameters.m_SignalGen.m_MaskImage); duplicator->Update(); m_Parameters.m_SignalGen.m_MaskImage = duplicator->GetOutput(); } if (m_Controls->m_TemplateComboBox->GetSelectedNode().IsNotNull() && mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( m_Controls->m_TemplateComboBox->GetSelectedNode())) // use parameters of selected DWI { mitk::Image::Pointer dwi = dynamic_cast(m_Controls->m_TemplateComboBox->GetSelectedNode()->GetData()); ItkDwiType::Pointer itkVectorImagePointer = ItkDwiType::New(); mitk::CastToItkImage(dwi, itkVectorImagePointer); m_Parameters.m_SignalGen.m_ImageRegion = itkVectorImagePointer->GetLargestPossibleRegion(); m_Parameters.m_SignalGen.m_ImageSpacing = itkVectorImagePointer->GetSpacing(); m_Parameters.m_SignalGen.m_ImageOrigin = itkVectorImagePointer->GetOrigin(); m_Parameters.m_SignalGen.m_ImageDirection = itkVectorImagePointer->GetDirection(); m_Parameters.SetBvalue(mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi)); m_Parameters.SetGradienDirections(mitk::DiffusionPropertyHelper::GetOriginalGradientContainer(dwi)); } else if (m_Controls->m_TemplateComboBox->GetSelectedNode().IsNotNull()) // use geometry of selected image { mitk::Image::Pointer img = dynamic_cast(m_Controls->m_TemplateComboBox->GetSelectedNode()->GetData()); itk::Image< float, 3 >::Pointer itkImg = itk::Image< float, 3 >::New(); CastToItkImage< itk::Image< float, 3 > >(img, itkImg); m_Parameters.m_SignalGen.m_ImageRegion = itkImg->GetLargestPossibleRegion(); m_Parameters.m_SignalGen.m_ImageSpacing = itkImg->GetSpacing(); m_Parameters.m_SignalGen.m_ImageOrigin = itkImg->GetOrigin(); m_Parameters.m_SignalGen.m_ImageDirection = itkImg->GetDirection(); if (m_Controls->m_UseBvalsBvecsBox->isChecked()) { double bval; m_Parameters.SetGradienDirections( mitk::gradients::ReadBvalsBvecs(m_Controls->m_LoadBvalsEdit->text().toStdString(), m_Controls->m_LoadBvecsEdit->text().toStdString(), bval) ); m_Parameters.SetBvalue(bval); } else { m_Parameters.SetNumWeightedVolumes(m_Controls->m_NumGradientsBox->value()); m_Parameters.SetBvalue(m_Controls->m_BvalueBox->value()); m_Parameters.GenerateGradientHalfShell(); } } else if (m_Parameters.m_SignalGen.m_MaskImage.IsNotNull()) // use geometry of mask image { ItkUcharImgType::Pointer itkImg = m_Parameters.m_SignalGen.m_MaskImage; m_Parameters.m_SignalGen.m_ImageRegion = itkImg->GetLargestPossibleRegion(); m_Parameters.m_SignalGen.m_ImageSpacing = itkImg->GetSpacing(); m_Parameters.m_SignalGen.m_ImageOrigin = itkImg->GetOrigin(); m_Parameters.m_SignalGen.m_ImageDirection = itkImg->GetDirection(); if (m_Controls->m_UseBvalsBvecsBox->isChecked()) { double bval; m_Parameters.SetGradienDirections( mitk::gradients::ReadBvalsBvecs(m_Controls->m_LoadBvalsEdit->text().toStdString(), m_Controls->m_LoadBvecsEdit->text().toStdString(), bval) ); m_Parameters.SetBvalue(bval); } else { m_Parameters.SetNumWeightedVolumes(m_Controls->m_NumGradientsBox->value()); m_Parameters.SetBvalue(m_Controls->m_BvalueBox->value()); m_Parameters.GenerateGradientHalfShell(); } } else // use GUI parameters { m_Parameters.m_SignalGen.m_ImageRegion.SetSize(0, m_Controls->m_SizeX->value()); m_Parameters.m_SignalGen.m_ImageRegion.SetSize(1, m_Controls->m_SizeY->value()); m_Parameters.m_SignalGen.m_ImageRegion.SetSize(2, m_Controls->m_SizeZ->value()); m_Parameters.m_SignalGen.m_ImageSpacing[0] = m_Controls->m_SpacingX->value(); m_Parameters.m_SignalGen.m_ImageSpacing[1] = m_Controls->m_SpacingY->value(); m_Parameters.m_SignalGen.m_ImageSpacing[2] = m_Controls->m_SpacingZ->value(); m_Parameters.m_SignalGen.m_ImageOrigin[0] = m_Parameters.m_SignalGen.m_ImageSpacing[0]/2; m_Parameters.m_SignalGen.m_ImageOrigin[1] = m_Parameters.m_SignalGen.m_ImageSpacing[1]/2; m_Parameters.m_SignalGen.m_ImageOrigin[2] = m_Parameters.m_SignalGen.m_ImageSpacing[2]/2; m_Parameters.m_SignalGen.m_ImageDirection.SetIdentity(); if (m_Controls->m_UseBvalsBvecsBox->isChecked()) { double bval; m_Parameters.SetGradienDirections( mitk::gradients::ReadBvalsBvecs(m_Controls->m_LoadBvalsEdit->text().toStdString(), m_Controls->m_LoadBvecsEdit->text().toStdString(), bval) ); m_Parameters.SetBvalue(bval); } else { m_Parameters.SetNumWeightedVolumes(m_Controls->m_NumGradientsBox->value()); m_Parameters.SetBvalue(m_Controls->m_BvalueBox->value()); m_Parameters.GenerateGradientHalfShell(); } } // signal relaxation m_Parameters.m_SignalGen.m_DoSimulateRelaxation = false; if (m_Controls->m_RelaxationBox->isChecked()) { m_Parameters.m_SignalGen.m_DoSimulateRelaxation = true; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Relaxation", BoolProperty::New(true)); m_Parameters.m_Misc.m_ArtifactModelString += "_RELAX"; } m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition = m_Parameters.m_SignalGen.m_DoSimulateRelaxation; // N/2 ghosts m_Parameters.m_Misc.m_DoAddGhosts = m_Controls->m_AddGhosts->isChecked(); m_Parameters.m_SignalGen.m_KspaceLineOffset = m_Controls->m_kOffsetBox->value(); if (m_Controls->m_AddGhosts->isChecked()) { m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; m_Parameters.m_Misc.m_ArtifactModelString += "_GHOST"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Ghost", DoubleProperty::New(m_Parameters.m_SignalGen.m_KspaceLineOffset)); } // Aliasing m_Parameters.m_Misc.m_DoAddAliasing = m_Controls->m_AddAliasing->isChecked(); m_Parameters.m_SignalGen.m_CroppingFactor = (100-m_Controls->m_WrapBox->value())/100; if (m_Controls->m_AddAliasing->isChecked()) { m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; m_Parameters.m_Misc.m_ArtifactModelString += "_ALIASING"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Aliasing", DoubleProperty::New(m_Controls->m_WrapBox->value())); } // Spikes m_Parameters.m_Misc.m_DoAddSpikes = m_Controls->m_AddSpikes->isChecked(); m_Parameters.m_SignalGen.m_Spikes = m_Controls->m_SpikeNumBox->value(); m_Parameters.m_SignalGen.m_SpikeAmplitude = m_Controls->m_SpikeScaleBox->value(); if (m_Controls->m_AddSpikes->isChecked()) { m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; m_Parameters.m_Misc.m_ArtifactModelString += "_SPIKES"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Spikes.Number", IntProperty::New(m_Parameters.m_SignalGen.m_Spikes)); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Spikes.Amplitude", DoubleProperty::New(m_Parameters.m_SignalGen.m_SpikeAmplitude)); } // Drift m_Parameters.m_SignalGen.m_DoAddDrift = m_Controls->m_AddDrift->isChecked(); - m_Parameters.m_SignalGen.m_Drift = m_Controls->m_DriftFactor->value()/100; + m_Parameters.m_SignalGen.m_Drift = static_cast(m_Controls->m_DriftFactor->value())/100; if (m_Controls->m_AddDrift->isChecked()) { m_Parameters.m_Misc.m_ArtifactModelString += "_DRIFT"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Drift", FloatProperty::New(m_Parameters.m_SignalGen.m_Drift)); } // gibbs ringing m_Parameters.m_SignalGen.m_DoAddGibbsRinging = m_Controls->m_AddGibbsRinging->isChecked(); + m_Parameters.m_SignalGen.m_ZeroRinging = m_Controls->m_ZeroRinging->value(); if (m_Controls->m_AddGibbsRinging->isChecked()) { m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Ringing", BoolProperty::New(true)); m_Parameters.m_Misc.m_ArtifactModelString += "_RINGING"; } // add distortions m_Parameters.m_Misc.m_DoAddDistortions = m_Controls->m_AddDistortions->isChecked(); if (m_Controls->m_AddDistortions->isChecked() && m_Controls->m_FrequencyMapBox->GetSelectedNode().IsNotNull()) { mitk::DataNode::Pointer fMapNode = m_Controls->m_FrequencyMapBox->GetSelectedNode(); mitk::Image* img = dynamic_cast(fMapNode->GetData()); ItkFloatImgType::Pointer itkImg = ItkFloatImgType::New(); CastToItkImage< ItkFloatImgType >(img, itkImg); if (m_Controls->m_TemplateComboBox->GetSelectedNode().IsNull()) // use geometry of frequency map { m_Parameters.m_SignalGen.m_ImageRegion = itkImg->GetLargestPossibleRegion(); m_Parameters.m_SignalGen.m_ImageSpacing = itkImg->GetSpacing(); m_Parameters.m_SignalGen.m_ImageOrigin = itkImg->GetOrigin(); m_Parameters.m_SignalGen.m_ImageDirection = itkImg->GetDirection(); } m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; itk::ImageDuplicator::Pointer duplicator = itk::ImageDuplicator::New(); duplicator->SetInputImage(itkImg); duplicator->Update(); m_Parameters.m_SignalGen.m_FrequencyMap = duplicator->GetOutput(); m_Parameters.m_Misc.m_ArtifactModelString += "_DISTORTED"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Distortions", BoolProperty::New(true)); } m_Parameters.m_SignalGen.m_EddyStrength = m_Controls->m_EddyGradientStrength->value(); m_Parameters.m_Misc.m_DoAddEddyCurrents = m_Controls->m_AddEddy->isChecked(); if (m_Controls->m_AddEddy->isChecked()) { m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; m_Parameters.m_Misc.m_ArtifactModelString += "_EDDY"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Eddy-strength", DoubleProperty::New(m_Parameters.m_SignalGen.m_EddyStrength)); } // Motion m_Parameters.m_SignalGen.m_DoAddMotion = false; m_Parameters.m_SignalGen.m_DoRandomizeMotion = m_Controls->m_RandomMotion->isChecked(); m_Parameters.m_SignalGen.m_Translation[0] = m_Controls->m_MaxTranslationBoxX->value(); m_Parameters.m_SignalGen.m_Translation[1] = m_Controls->m_MaxTranslationBoxY->value(); m_Parameters.m_SignalGen.m_Translation[2] = m_Controls->m_MaxTranslationBoxZ->value(); m_Parameters.m_SignalGen.m_Rotation[0] = m_Controls->m_MaxRotationBoxX->value(); m_Parameters.m_SignalGen.m_Rotation[1] = m_Controls->m_MaxRotationBoxY->value(); m_Parameters.m_SignalGen.m_Rotation[2] = m_Controls->m_MaxRotationBoxZ->value(); m_Parameters.m_SignalGen.m_MotionVolumes.clear(); m_Parameters.m_Misc.m_MotionVolumesBox = m_Controls->m_MotionVolumesBox->text().toStdString(); if ( m_Controls->m_AddMotion->isChecked()) { m_Parameters.m_SignalGen.m_DoAddMotion = true; m_Parameters.m_Misc.m_ArtifactModelString += "_MOTION"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Random", BoolProperty::New(m_Parameters.m_SignalGen.m_DoRandomizeMotion)); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Translation-x", DoubleProperty::New(m_Parameters.m_SignalGen.m_Translation[0])); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Translation-y", DoubleProperty::New(m_Parameters.m_SignalGen.m_Translation[1])); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Translation-z", DoubleProperty::New(m_Parameters.m_SignalGen.m_Translation[2])); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Rotation-x", DoubleProperty::New(m_Parameters.m_SignalGen.m_Rotation[0])); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Rotation-y", DoubleProperty::New(m_Parameters.m_SignalGen.m_Rotation[1])); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Rotation-z", DoubleProperty::New(m_Parameters.m_SignalGen.m_Rotation[2])); if ( m_Parameters.m_Misc.m_MotionVolumesBox == "random" ) { for ( size_t i=0; i < m_Parameters.m_SignalGen.GetNumVolumes(); ++i ) { m_Parameters.m_SignalGen.m_MotionVolumes.push_back( bool( rand()%2 ) ); } MITK_DEBUG << "QmitkFiberfoxView.cpp: Case m_Misc.m_MotionVolumesBox == \"random\"."; } else if ( ! m_Parameters.m_Misc.m_MotionVolumesBox.empty() ) { std::stringstream stream( m_Parameters.m_Misc.m_MotionVolumesBox ); std::vector numbers; int number = std::numeric_limits::max(); while( stream >> number ) { if( number < std::numeric_limits::max() ) { numbers.push_back( number ); } } // If a list of negative numbers is given: if( *(std::min_element( numbers.begin(), numbers.end() )) < 0 && *(std::max_element( numbers.begin(), numbers.end() )) <= 0 ) // cave: -0 == +0 { for ( size_t i=0; i < m_Parameters.m_SignalGen.GetNumVolumes(); ++i ) { m_Parameters.m_SignalGen.m_MotionVolumes.push_back( true ); } // set all true except those given. for( auto iter = std::begin( numbers ); iter != std::end( numbers ); ++iter ) { if ( -(*iter) < (int)m_Parameters.m_SignalGen.GetNumVolumes() && -(*iter) >= 0 ) { m_Parameters.m_SignalGen.m_MotionVolumes.at( -(*iter) ) = false; } } MITK_DEBUG << "QmitkFiberfoxView.cpp: Case list of negative numbers."; } // If a list of positive numbers is given: else if( *(std::min_element( numbers.begin(), numbers.end() )) >= 0 && *(std::max_element( numbers.begin(), numbers.end() )) >= 0 ) { for ( size_t i=0; i < m_Parameters.m_SignalGen.GetNumVolumes(); ++i ) { m_Parameters.m_SignalGen.m_MotionVolumes.push_back( false ); } // set all false except those given. for( auto iter = std::begin( numbers ); iter != std::end( numbers ); ++iter ) { if ( *iter < (int)m_Parameters.m_SignalGen.GetNumVolumes() && *iter >= 0 ) { m_Parameters.m_SignalGen.m_MotionVolumes.at( *iter ) = true; } } MITK_DEBUG << "QmitkFiberfoxView.cpp: Case list of positive numbers."; } else { MITK_ERROR << "QmitkFiberfoxView.cpp: Inconsistent list of numbers in m_MotionVolumesBox."; } } else { m_Parameters.m_Misc.m_MotionVolumesBox = ""; // set empty. m_Controls->m_MotionVolumesBox->setText(""); for (unsigned int i=0; im_AcquisitionTypeBox->currentIndex(); m_Parameters.m_SignalGen.m_CoilSensitivityProfile = (SignalGenerationParameters::CoilSensitivityProfile)m_Controls->m_CoilSensBox->currentIndex(); m_Parameters.m_SignalGen.m_NumberOfCoils = m_Controls->m_NumCoilsBox->value(); m_Parameters.m_SignalGen.m_PartialFourier = m_Controls->m_PartialFourier->value(); m_Parameters.m_SignalGen.m_ReversePhase = m_Controls->m_ReversePhaseBox->isChecked(); m_Parameters.m_SignalGen.m_tLine = m_Controls->m_LineReadoutTimeBox->value(); m_Parameters.m_SignalGen.m_tInhom = m_Controls->m_T2starBox->value(); m_Parameters.m_SignalGen.m_tEcho = m_Controls->m_TEbox->value(); m_Parameters.m_SignalGen.m_tRep = m_Controls->m_TRbox->value(); m_Parameters.m_SignalGen.m_DoDisablePartialVolume = m_Controls->m_EnforcePureFiberVoxelsBox->isChecked(); m_Parameters.m_SignalGen.m_AxonRadius = m_Controls->m_FiberRadius->value(); m_Parameters.m_SignalGen.m_SignalScale = m_Controls->m_SignalScaleBox->value(); double voxelVolume = m_Parameters.m_SignalGen.m_ImageSpacing[0] * m_Parameters.m_SignalGen.m_ImageSpacing[1] * m_Parameters.m_SignalGen.m_ImageSpacing[2]; if ( m_Parameters.m_SignalGen.m_SignalScale*voxelVolume > itk::NumericTraits::max()*0.75 ) { m_Parameters.m_SignalGen.m_SignalScale = itk::NumericTraits::max()*0.75/voxelVolume; m_Controls->m_SignalScaleBox->setValue(m_Parameters.m_SignalGen.m_SignalScale); QMessageBox::information( nullptr, "Warning", "Maximum signal exceeding data type limits. Automatically adjusted to " + QString::number(m_Parameters.m_SignalGen.m_SignalScale) + " to obtain a maximum signal of 75% of the data type maximum." " Relaxation and other effects that affect the signal intensities are not accounted for."); } // Noise m_Parameters.m_Misc.m_DoAddNoise = m_Controls->m_AddNoise->isChecked(); m_Parameters.m_SignalGen.m_NoiseVariance = m_Controls->m_NoiseLevel->value(); if (m_Controls->m_AddNoise->isChecked()) { switch (m_Controls->m_NoiseDistributionBox->currentIndex()) { case 0: { if (m_Parameters.m_SignalGen.m_NoiseVariance>0) { m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; m_Parameters.m_Misc.m_ArtifactModelString += "_COMPLEX-GAUSSIAN-"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Noise-Distribution", StringProperty::New("Complex Gaussian")); } break; } case 1: { if (m_Parameters.m_SignalGen.m_NoiseVariance>0) { m_Parameters.m_NoiseModel = std::make_shared< mitk::RicianNoiseModel >(); m_Parameters.m_Misc.m_ArtifactModelString += "_RICIAN-"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Noise-Distribution", StringProperty::New("Rician")); m_Parameters.m_NoiseModel->SetNoiseVariance(m_Parameters.m_SignalGen.m_NoiseVariance); } break; } case 2: { if (m_Parameters.m_SignalGen.m_NoiseVariance>0) { m_Parameters.m_NoiseModel = std::make_shared< mitk::ChiSquareNoiseModel >(); m_Parameters.m_Misc.m_ArtifactModelString += "_CHISQUARED-"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Noise-Distribution", StringProperty::New("Chi-squared")); m_Parameters.m_NoiseModel->SetNoiseVariance(m_Parameters.m_SignalGen.m_NoiseVariance); } break; } default: { if (m_Parameters.m_SignalGen.m_NoiseVariance>0) { m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; m_Parameters.m_Misc.m_ArtifactModelString += "_COMPLEX-GAUSSIAN-"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Noise-Distribution", StringProperty::New("Complex Gaussian")); } break; } } if (m_Parameters.m_SignalGen.m_NoiseVariance>0) { m_Parameters.m_Misc.m_ArtifactModelString += QString::number(m_Parameters.m_SignalGen.m_NoiseVariance).toStdString(); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Noise-Variance", DoubleProperty::New(m_Parameters.m_SignalGen.m_NoiseVariance)); } } // signal models { // compartment 1 switch (m_Controls->m_Compartment1Box->currentIndex()) { case 0: { mitk::StickModel* model = new mitk::StickModel(); model->SetGradientList(m_Parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(m_Parameters.m_SignalGen.GetBvalue()); model->SetDiffusivity(m_Controls->m_StickWidget1->GetD()); model->SetT2(m_Controls->m_StickWidget1->GetT2()); model->SetT1(m_Controls->m_StickWidget1->GetT1()); model->m_CompartmentId = 1; m_Parameters.m_FiberModelList.push_back(model); m_Parameters.m_Misc.m_SignalModelString += "Stick"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Stick") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D", DoubleProperty::New(m_Controls->m_StickWidget1->GetD()) ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.T2", DoubleProperty::New(model->GetT2()) ); break; } case 1: { mitk::TensorModel* model = new mitk::TensorModel(); model->SetGradientList(m_Parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(m_Parameters.m_SignalGen.GetBvalue()); model->SetDiffusivity1(m_Controls->m_ZeppelinWidget1->GetD1()); model->SetDiffusivity2(m_Controls->m_ZeppelinWidget1->GetD2()); model->SetDiffusivity3(m_Controls->m_ZeppelinWidget1->GetD2()); model->SetT2(m_Controls->m_ZeppelinWidget1->GetT2()); model->SetT1(m_Controls->m_ZeppelinWidget1->GetT1()); model->m_CompartmentId = 1; m_Parameters.m_FiberModelList.push_back(model); m_Parameters.m_Misc.m_SignalModelString += "Zeppelin"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Zeppelin") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D1", DoubleProperty::New(m_Controls->m_ZeppelinWidget1->GetD1()) ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D2", DoubleProperty::New(m_Controls->m_ZeppelinWidget1->GetD2()) ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.T2", DoubleProperty::New(model->GetT2()) ); break; } case 2: { mitk::TensorModel* model = new mitk::TensorModel(); model->SetGradientList(m_Parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(m_Parameters.m_SignalGen.GetBvalue()); model->SetDiffusivity1(m_Controls->m_TensorWidget1->GetD1()); model->SetDiffusivity2(m_Controls->m_TensorWidget1->GetD2()); model->SetDiffusivity3(m_Controls->m_TensorWidget1->GetD3()); model->SetT2(m_Controls->m_TensorWidget1->GetT2()); model->SetT1(m_Controls->m_TensorWidget1->GetT1()); model->m_CompartmentId = 1; m_Parameters.m_FiberModelList.push_back(model); m_Parameters.m_Misc.m_SignalModelString += "Tensor"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Tensor") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D1", DoubleProperty::New(m_Controls->m_TensorWidget1->GetD1()) ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D2", DoubleProperty::New(m_Controls->m_TensorWidget1->GetD2()) ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D3", DoubleProperty::New(m_Controls->m_TensorWidget1->GetD3()) ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.T2", DoubleProperty::New(model->GetT2()) ); break; } case 3: { mitk::RawShModel* model = new mitk::RawShModel(); m_Parameters.m_SignalGen.m_DoSimulateRelaxation = false; model->SetGradientList(m_Parameters.m_SignalGen.GetGradientDirections()); model->SetMaxNumKernels(m_Controls->m_PrototypeWidget1->GetNumberOfSamples()); model->SetFaRange(m_Controls->m_PrototypeWidget1->GetMinFa(), m_Controls->m_PrototypeWidget1->GetMaxFa()); model->SetAdcRange(m_Controls->m_PrototypeWidget1->GetMinAdc(), m_Controls->m_PrototypeWidget1->GetMaxAdc()); model->m_CompartmentId = 1; m_Parameters.m_FiberModelList.push_back(model); m_Parameters.m_Misc.m_SignalModelString += "Prototype"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Prototype") ); break; } } if (m_Controls->m_Comp1VolumeFraction->GetSelectedNode().IsNotNull()) { mitk::DataNode::Pointer volumeNode = m_Controls->m_Comp1VolumeFraction->GetSelectedNode(); ItkDoubleImgType::Pointer comp1VolumeImage = ItkDoubleImgType::New(); mitk::Image* img = dynamic_cast(volumeNode->GetData()); CastToItkImage< ItkDoubleImgType >(img, comp1VolumeImage); m_Parameters.m_FiberModelList.back()->SetVolumeFractionImage(comp1VolumeImage); } // compartment 2 switch (m_Controls->m_Compartment2Box->currentIndex()) { case 0: break; case 1: { mitk::StickModel* model = new mitk::StickModel(); model->SetGradientList(m_Parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(m_Parameters.m_SignalGen.GetBvalue()); model->SetDiffusivity(m_Controls->m_StickWidget2->GetD()); model->SetT2(m_Controls->m_StickWidget2->GetT2()); model->SetT1(m_Controls->m_StickWidget2->GetT1()); model->m_CompartmentId = 2; m_Parameters.m_FiberModelList.push_back(model); m_Parameters.m_Misc.m_SignalModelString += "Stick"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Description", StringProperty::New("Inter-axonal compartment") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Model", StringProperty::New("Stick") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D", DoubleProperty::New(m_Controls->m_StickWidget2->GetD()) ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.T2", DoubleProperty::New(model->GetT2()) ); break; } case 2: { mitk::TensorModel* model = new mitk::TensorModel(); model->SetGradientList(m_Parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(m_Parameters.m_SignalGen.GetBvalue()); model->SetDiffusivity1(m_Controls->m_ZeppelinWidget2->GetD1()); model->SetDiffusivity2(m_Controls->m_ZeppelinWidget2->GetD2()); model->SetDiffusivity3(m_Controls->m_ZeppelinWidget2->GetD2()); model->SetT2(m_Controls->m_ZeppelinWidget2->GetT2()); model->SetT1(m_Controls->m_ZeppelinWidget2->GetT1()); model->m_CompartmentId = 2; m_Parameters.m_FiberModelList.push_back(model); m_Parameters.m_Misc.m_SignalModelString += "Zeppelin"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Description", StringProperty::New("Inter-axonal compartment") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Model", StringProperty::New("Zeppelin") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D1", DoubleProperty::New(m_Controls->m_ZeppelinWidget2->GetD1()) ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D2", DoubleProperty::New(m_Controls->m_ZeppelinWidget2->GetD2()) ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.T2", DoubleProperty::New(model->GetT2()) ); break; } case 3: { mitk::TensorModel* model = new mitk::TensorModel(); model->SetGradientList(m_Parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(m_Parameters.m_SignalGen.GetBvalue()); model->SetDiffusivity1(m_Controls->m_TensorWidget2->GetD1()); model->SetDiffusivity2(m_Controls->m_TensorWidget2->GetD2()); model->SetDiffusivity3(m_Controls->m_TensorWidget2->GetD3()); model->SetT2(m_Controls->m_TensorWidget2->GetT2()); model->SetT1(m_Controls->m_TensorWidget2->GetT1()); model->m_CompartmentId = 2; m_Parameters.m_FiberModelList.push_back(model); m_Parameters.m_Misc.m_SignalModelString += "Tensor"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Description", StringProperty::New("Inter-axonal compartment") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Model", StringProperty::New("Tensor") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D1", DoubleProperty::New(m_Controls->m_TensorWidget2->GetD1()) ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D2", DoubleProperty::New(m_Controls->m_TensorWidget2->GetD2()) ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D3", DoubleProperty::New(m_Controls->m_TensorWidget2->GetD3()) ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.T2", DoubleProperty::New(model->GetT2()) ); break; } } if (m_Controls->m_Comp2VolumeFraction->GetSelectedNode().IsNotNull() && m_Parameters.m_FiberModelList.size()==2) { mitk::DataNode::Pointer volumeNode = m_Controls->m_Comp2VolumeFraction->GetSelectedNode(); ItkDoubleImgType::Pointer comp1VolumeImage = ItkDoubleImgType::New(); mitk::Image* img = dynamic_cast(volumeNode->GetData()); CastToItkImage< ItkDoubleImgType >(img, comp1VolumeImage); m_Parameters.m_FiberModelList.back()->SetVolumeFractionImage(comp1VolumeImage); } // compartment 3 switch (m_Controls->m_Compartment3Box->currentIndex()) { case 0: { mitk::BallModel* model = new mitk::BallModel(); model->SetGradientList(m_Parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(m_Parameters.m_SignalGen.GetBvalue()); model->SetDiffusivity(m_Controls->m_BallWidget1->GetD()); model->SetT2(m_Controls->m_BallWidget1->GetT2()); model->SetT1(m_Controls->m_BallWidget1->GetT1()); model->m_CompartmentId = 3; m_Parameters.m_NonFiberModelList.push_back(model); m_Parameters.m_Misc.m_SignalModelString += "Ball"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Ball") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.D", DoubleProperty::New(m_Controls->m_BallWidget1->GetD()) ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.T2", DoubleProperty::New(model->GetT2()) ); break; } case 1: { mitk::AstroStickModel* model = new mitk::AstroStickModel(); model->SetGradientList(m_Parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(m_Parameters.m_SignalGen.GetBvalue()); model->SetDiffusivity(m_Controls->m_AstrosticksWidget1->GetD()); model->SetT2(m_Controls->m_AstrosticksWidget1->GetT2()); model->SetT1(m_Controls->m_AstrosticksWidget1->GetT1()); model->SetRandomizeSticks(m_Controls->m_AstrosticksWidget1->GetRandomizeSticks()); model->m_CompartmentId = 3; m_Parameters.m_NonFiberModelList.push_back(model); m_Parameters.m_Misc.m_SignalModelString += "Astrosticks"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Astrosticks") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.D", DoubleProperty::New(m_Controls->m_AstrosticksWidget1->GetD()) ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.T2", DoubleProperty::New(model->GetT2()) ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.RandomSticks", BoolProperty::New(m_Controls->m_AstrosticksWidget1->GetRandomizeSticks()) ); break; } case 2: { mitk::DotModel* model = new mitk::DotModel(); model->SetGradientList(m_Parameters.m_SignalGen.GetGradientDirections()); model->SetT2(m_Controls->m_DotWidget1->GetT2()); model->SetT1(m_Controls->m_DotWidget1->GetT1()); model->m_CompartmentId = 3; m_Parameters.m_NonFiberModelList.push_back(model); m_Parameters.m_Misc.m_SignalModelString += "Dot"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Dot") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.T2", DoubleProperty::New(model->GetT2()) ); break; } case 3: { mitk::RawShModel* model = new mitk::RawShModel(); m_Parameters.m_SignalGen.m_DoSimulateRelaxation = false; model->SetGradientList(m_Parameters.m_SignalGen.GetGradientDirections()); model->SetMaxNumKernels(m_Controls->m_PrototypeWidget3->GetNumberOfSamples()); model->SetFaRange(m_Controls->m_PrototypeWidget3->GetMinFa(), m_Controls->m_PrototypeWidget3->GetMaxFa()); model->SetAdcRange(m_Controls->m_PrototypeWidget3->GetMinAdc(), m_Controls->m_PrototypeWidget3->GetMaxAdc()); model->m_CompartmentId = 3; m_Parameters.m_NonFiberModelList.push_back(model); m_Parameters.m_Misc.m_SignalModelString += "Prototype"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Prototype") ); break; } } if (m_Controls->m_Comp3VolumeFraction->GetSelectedNode().IsNotNull()) { mitk::DataNode::Pointer volumeNode = m_Controls->m_Comp3VolumeFraction->GetSelectedNode(); ItkDoubleImgType::Pointer comp1VolumeImage = ItkDoubleImgType::New(); mitk::Image* img = dynamic_cast(volumeNode->GetData()); CastToItkImage< ItkDoubleImgType >(img, comp1VolumeImage); m_Parameters.m_NonFiberModelList.back()->SetVolumeFractionImage(comp1VolumeImage); } switch (m_Controls->m_Compartment4Box->currentIndex()) { case 0: break; case 1: { mitk::BallModel* model = new mitk::BallModel(); model->SetGradientList(m_Parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(m_Parameters.m_SignalGen.GetBvalue()); model->SetDiffusivity(m_Controls->m_BallWidget2->GetD()); model->SetT2(m_Controls->m_BallWidget2->GetT2()); model->SetT1(m_Controls->m_BallWidget2->GetT1()); model->m_CompartmentId = 4; m_Parameters.m_NonFiberModelList.push_back(model); m_Parameters.m_Misc.m_SignalModelString += "Ball"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Ball") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.D", DoubleProperty::New(m_Controls->m_BallWidget2->GetD()) ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.T2", DoubleProperty::New(model->GetT2()) ); break; } case 2: { mitk::AstroStickModel* model = new mitk::AstroStickModel(); model->SetGradientList(m_Parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(m_Parameters.m_SignalGen.GetBvalue()); model->SetDiffusivity(m_Controls->m_AstrosticksWidget2->GetD()); model->SetT2(m_Controls->m_AstrosticksWidget2->GetT2()); model->SetT1(m_Controls->m_AstrosticksWidget2->GetT1()); model->SetRandomizeSticks(m_Controls->m_AstrosticksWidget2->GetRandomizeSticks()); model->m_CompartmentId = 4; m_Parameters.m_NonFiberModelList.push_back(model); m_Parameters.m_Misc.m_SignalModelString += "Astrosticks"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Astrosticks") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.D", DoubleProperty::New(m_Controls->m_AstrosticksWidget2->GetD()) ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.T2", DoubleProperty::New(model->GetT2()) ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.RandomSticks", BoolProperty::New(m_Controls->m_AstrosticksWidget2->GetRandomizeSticks()) ); break; } case 3: { mitk::DotModel* model = new mitk::DotModel(); model->SetGradientList(m_Parameters.m_SignalGen.GetGradientDirections()); model->SetT2(m_Controls->m_DotWidget2->GetT2()); model->SetT1(m_Controls->m_DotWidget2->GetT1()); model->m_CompartmentId = 4; m_Parameters.m_NonFiberModelList.push_back(model); m_Parameters.m_Misc.m_SignalModelString += "Dot"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Dot") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.T2", DoubleProperty::New(model->GetT2()) ); break; } case 4: { mitk::RawShModel* model = new mitk::RawShModel(); m_Parameters.m_SignalGen.m_DoSimulateRelaxation = false; model->SetGradientList(m_Parameters.m_SignalGen.GetGradientDirections()); model->SetMaxNumKernels(m_Controls->m_PrototypeWidget4->GetNumberOfSamples()); model->SetFaRange(m_Controls->m_PrototypeWidget4->GetMinFa(), m_Controls->m_PrototypeWidget4->GetMaxFa()); model->SetAdcRange(m_Controls->m_PrototypeWidget4->GetMinAdc(), m_Controls->m_PrototypeWidget4->GetMaxAdc()); model->m_CompartmentId = 4; m_Parameters.m_NonFiberModelList.push_back(model); m_Parameters.m_Misc.m_SignalModelString += "Prototype"; m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Prototype") ); break; } } if (m_Controls->m_Comp4VolumeFraction->GetSelectedNode().IsNotNull() && m_Parameters.m_NonFiberModelList.size()==2) { mitk::DataNode::Pointer volumeNode = m_Controls->m_Comp4VolumeFraction->GetSelectedNode(); ItkDoubleImgType::Pointer compVolumeImage = ItkDoubleImgType::New(); mitk::Image* img = dynamic_cast(volumeNode->GetData()); CastToItkImage< ItkDoubleImgType >(img, compVolumeImage); m_Parameters.m_NonFiberModelList.back()->SetVolumeFractionImage(compVolumeImage); } } m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.SignalScale", IntProperty::New(m_Parameters.m_SignalGen.m_SignalScale)); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.FiberRadius", IntProperty::New(m_Parameters.m_SignalGen.m_AxonRadius)); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Tinhom", DoubleProperty::New(m_Parameters.m_SignalGen.m_tInhom)); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Tline", DoubleProperty::New(m_Parameters.m_SignalGen.m_tLine)); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.TE", DoubleProperty::New(m_Parameters.m_SignalGen.m_tEcho)); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.b-value", DoubleProperty::New(m_Parameters.m_SignalGen.GetBvalue())); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.NoPartialVolume", BoolProperty::New(m_Parameters.m_SignalGen.m_DoDisablePartialVolume)); m_Parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Relaxation", BoolProperty::New(m_Parameters.m_SignalGen.m_DoSimulateRelaxation)); m_Parameters.m_Misc.m_ResultNode->AddProperty("binary", BoolProperty::New(false)); } void QmitkFiberfoxView::SaveParameters(QString filename) { UpdateParametersFromGui(); std::vector< int > bVals = m_Parameters.m_SignalGen.GetBvalues(); std::cout << "b-values: "; for (auto v : bVals) std::cout << v << " "; std::cout << std::endl; bool ok = true; bool first = true; bool dosampling = false; mitk::Image::Pointer diffImg = nullptr; itk::Image< itk::DiffusionTensor3D< double >, 3 >::Pointer tensorImage = nullptr; const int shOrder = 2; typedef itk::AnalyticalDiffusionQballReconstructionImageFilter QballFilterType; QballFilterType::CoefficientImageType::Pointer itkFeatureImage = nullptr; ItkDoubleImgType::Pointer adcImage = nullptr; for (unsigned int i=0; i* model = nullptr; if (i* >(m_Parameters.m_FiberModelList.at(i)); } else { model = dynamic_cast< mitk::RawShModel<>* >(m_Parameters.m_NonFiberModelList.at(i-m_Parameters.m_FiberModelList.size())); } if ( model!=nullptr && model->GetNumberOfKernels() <= 0 ) { if (first==true) { if ( QMessageBox::question(nullptr, "Prototype signal sampling", "Do you want to sample prototype signals from the selected diffusion-weighted imag and save them?", QMessageBox::Yes, QMessageBox::No) == QMessageBox::Yes ) dosampling = true; first = false; if ( dosampling && (m_Controls->m_TemplateComboBox->GetSelectedNode().IsNull() || !mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_Controls->m_TemplateComboBox->GetSelectedNode()->GetData()) ) ) ) { QMessageBox::information(nullptr, "Parameter file not saved", "No diffusion-weighted image selected to sample signal from."); return; } else if (dosampling) { diffImg = dynamic_cast(m_Controls->m_TemplateComboBox->GetSelectedNode()->GetData()); typedef itk::DiffusionTensor3DReconstructionImageFilter< short, short, double > TensorReconstructionImageFilterType; TensorReconstructionImageFilterType::Pointer filter = TensorReconstructionImageFilterType::New(); ItkDwiType::Pointer itkVectorImagePointer = ItkDwiType::New(); mitk::CastToItkImage(diffImg, itkVectorImagePointer); filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(diffImg)); filter->SetGradientImage(mitk::DiffusionPropertyHelper::GetGradientContainer(diffImg), itkVectorImagePointer ); filter->Update(); tensorImage = filter->GetOutput(); QballFilterType::Pointer qballfilter = QballFilterType::New(); qballfilter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(diffImg)); qballfilter->SetGradientImage(mitk::DiffusionPropertyHelper::GetGradientContainer(diffImg), itkVectorImagePointer ); qballfilter->SetLambda(0.006); qballfilter->SetNormalizationMethod(QballFilterType::QBAR_RAW_SIGNAL); qballfilter->Update(); itkFeatureImage = qballfilter->GetCoefficientImage(); itk::AdcImageFilter< short, double >::Pointer adcFilter = itk::AdcImageFilter< short, double >::New(); adcFilter->SetInput( itkVectorImagePointer ); adcFilter->SetGradientDirections(mitk::DiffusionPropertyHelper::GetGradientContainer(diffImg)); adcFilter->SetB_value(mitk::DiffusionPropertyHelper::GetReferenceBValue(diffImg)); adcFilter->Update(); adcImage = adcFilter->GetOutput(); } } typedef itk::DiffusionTensor3DReconstructionImageFilter< short, short, double > TensorReconstructionImageFilterType; TensorReconstructionImageFilterType::Pointer filter = TensorReconstructionImageFilterType::New(); ItkDwiType::Pointer itkVectorImagePointer = ItkDwiType::New(); mitk::CastToItkImage(diffImg, itkVectorImagePointer); filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(diffImg)); filter->SetGradientImage(mitk::DiffusionPropertyHelper::GetGradientContainer(diffImg), itkVectorImagePointer ); filter->Update(); tensorImage = filter->GetOutput(); QballFilterType::Pointer qballfilter = QballFilterType::New(); qballfilter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(diffImg)); qballfilter->SetGradientImage(mitk::DiffusionPropertyHelper::GetGradientContainer(diffImg), itkVectorImagePointer ); qballfilter->SetLambda(0.006); qballfilter->SetNormalizationMethod(QballFilterType::QBAR_RAW_SIGNAL); qballfilter->Update(); itkFeatureImage = qballfilter->GetCoefficientImage(); itk::AdcImageFilter< short, double >::Pointer adcFilter = itk::AdcImageFilter< short, double >::New(); adcFilter->SetInput( itkVectorImagePointer ); adcFilter->SetGradientDirections(mitk::DiffusionPropertyHelper::GetGradientContainer(diffImg)); adcFilter->SetB_value(mitk::DiffusionPropertyHelper::GetReferenceBValue(diffImg)); adcFilter->Update(); adcImage = adcFilter->GetOutput(); if (dosampling && diffImg.IsNotNull()) { ok = model->SampleKernels(diffImg, m_Parameters.m_SignalGen.m_MaskImage, tensorImage, itkFeatureImage, adcImage); if (!ok) { QMessageBox::information( nullptr, "Parameter file not saved", "No valid prototype signals could be sampled."); return; } } } } m_Parameters.SaveParameters(filename.toStdString()); m_ParameterFile = filename; } void QmitkFiberfoxView::SaveParameters() { QString filename = QFileDialog::getSaveFileName( 0, tr("Save Parameters"), m_ParameterFile, tr("Fiberfox Parameters (*.ffp)") ); SaveParameters(filename); } void QmitkFiberfoxView::LoadParameters() { QString filename = QFileDialog::getOpenFileName(0, tr("Load Parameters"), QString(itksys::SystemTools::GetFilenamePath(m_ParameterFile.toStdString()).c_str()), tr("Fiberfox Parameters (*.ffp)") ); if(filename.isEmpty() || filename.isNull()) return; m_ParameterFile = filename; m_Parameters.LoadParameters(filename.toStdString()); if (m_Parameters.m_MissingTags.size()>0) { QString missing("Parameter file might be corrupted. The following parameters could not be read: "); missing += QString(m_Parameters.m_MissingTags.c_str()); missing += "\nDefault values have been assigned to the missing parameters."; QMessageBox::information( nullptr, "Warning!", missing); } // image generation parameters m_Controls->m_SizeX->setValue(m_Parameters.m_SignalGen.m_ImageRegion.GetSize(0)); m_Controls->m_SizeY->setValue(m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1)); m_Controls->m_SizeZ->setValue(m_Parameters.m_SignalGen.m_ImageRegion.GetSize(2)); m_Controls->m_SpacingX->setValue(m_Parameters.m_SignalGen.m_ImageSpacing[0]); m_Controls->m_SpacingY->setValue(m_Parameters.m_SignalGen.m_ImageSpacing[1]); m_Controls->m_SpacingZ->setValue(m_Parameters.m_SignalGen.m_ImageSpacing[2]); m_Controls->m_NumGradientsBox->setValue(m_Parameters.m_SignalGen.GetNumWeightedVolumes()); m_Controls->m_BvalueBox->setValue(m_Parameters.m_SignalGen.GetBvalue()); m_Controls->m_SignalScaleBox->setValue(m_Parameters.m_SignalGen.m_SignalScale); m_Controls->m_TEbox->setValue(m_Parameters.m_SignalGen.m_tEcho); m_Controls->m_LineReadoutTimeBox->setValue(m_Parameters.m_SignalGen.m_tLine); m_Controls->m_T2starBox->setValue(m_Parameters.m_SignalGen.m_tInhom); m_Controls->m_FiberRadius->setValue(m_Parameters.m_SignalGen.m_AxonRadius); m_Controls->m_RelaxationBox->setChecked(m_Parameters.m_SignalGen.m_DoSimulateRelaxation); m_Controls->m_EnforcePureFiberVoxelsBox->setChecked(m_Parameters.m_SignalGen.m_DoDisablePartialVolume); m_Controls->m_ReversePhaseBox->setChecked(m_Parameters.m_SignalGen.m_ReversePhase); m_Controls->m_PartialFourier->setValue(m_Parameters.m_SignalGen.m_PartialFourier); m_Controls->m_TRbox->setValue(m_Parameters.m_SignalGen.m_tRep); m_Controls->m_NumCoilsBox->setValue(m_Parameters.m_SignalGen.m_NumberOfCoils); m_Controls->m_CoilSensBox->setCurrentIndex(m_Parameters.m_SignalGen.m_CoilSensitivityProfile); m_Controls->m_AcquisitionTypeBox->setCurrentIndex(m_Parameters.m_SignalGen.m_AcquisitionType); if (!m_Parameters.m_Misc.m_BvalsFile.empty()) { m_Controls->m_UseBvalsBvecsBox->setChecked(true); m_Controls->m_LoadBvalsEdit->setText(QString(m_Parameters.m_Misc.m_BvalsFile.c_str())); } else m_Controls->m_LoadBvalsEdit->setText("-"); if (!m_Parameters.m_Misc.m_BvecsFile.empty()) { m_Controls->m_UseBvalsBvecsBox->setChecked(true); m_Controls->m_LoadBvecsEdit->setText(QString(m_Parameters.m_Misc.m_BvecsFile.c_str())); } else m_Controls->m_LoadBvecsEdit->setText("-"); if (m_Parameters.m_NoiseModel!=nullptr) { m_Controls->m_AddNoise->setChecked(m_Parameters.m_Misc.m_DoAddNoise); if (dynamic_cast*>(m_Parameters.m_NoiseModel.get())) { m_Controls->m_NoiseDistributionBox->setCurrentIndex(0); } else if (dynamic_cast*>(m_Parameters.m_NoiseModel.get())) { m_Controls->m_NoiseDistributionBox->setCurrentIndex(1); } m_Controls->m_NoiseLevel->setValue(m_Parameters.m_NoiseModel->GetNoiseVariance()); } else { m_Controls->m_AddNoise->setChecked(m_Parameters.m_Misc.m_DoAddNoise); m_Controls->m_NoiseLevel->setValue(m_Parameters.m_SignalGen.m_NoiseVariance); } m_Controls->m_VolumeFractionsBox->setChecked(m_Parameters.m_Misc.m_CheckOutputVolumeFractionsBox); m_Controls->m_AdvancedOptionsBox_2->setChecked(m_Parameters.m_Misc.m_CheckAdvancedSignalOptionsBox); m_Controls->m_AddGhosts->setChecked(m_Parameters.m_Misc.m_DoAddGhosts); m_Controls->m_AddAliasing->setChecked(m_Parameters.m_Misc.m_DoAddAliasing); m_Controls->m_AddDistortions->setChecked(m_Parameters.m_Misc.m_DoAddDistortions); m_Controls->m_AddSpikes->setChecked(m_Parameters.m_Misc.m_DoAddSpikes); m_Controls->m_AddEddy->setChecked(m_Parameters.m_Misc.m_DoAddEddyCurrents); m_Controls->m_AddDrift->setChecked(m_Parameters.m_SignalGen.m_DoAddDrift); m_Controls->m_kOffsetBox->setValue(m_Parameters.m_SignalGen.m_KspaceLineOffset); m_Controls->m_WrapBox->setValue(100*(1-m_Parameters.m_SignalGen.m_CroppingFactor)); m_Controls->m_DriftFactor->setValue(100*m_Parameters.m_SignalGen.m_Drift); m_Controls->m_SpikeNumBox->setValue(m_Parameters.m_SignalGen.m_Spikes); m_Controls->m_SpikeScaleBox->setValue(m_Parameters.m_SignalGen.m_SpikeAmplitude); m_Controls->m_EddyGradientStrength->setValue(m_Parameters.m_SignalGen.m_EddyStrength); m_Controls->m_AddGibbsRinging->setChecked(m_Parameters.m_SignalGen.m_DoAddGibbsRinging); + m_Controls->m_ZeroRinging->setValue(m_Parameters.m_SignalGen.m_ZeroRinging); m_Controls->m_AddMotion->setChecked(m_Parameters.m_SignalGen.m_DoAddMotion); m_Controls->m_RandomMotion->setChecked(m_Parameters.m_SignalGen.m_DoRandomizeMotion); m_Controls->m_MotionVolumesBox->setText(QString(m_Parameters.m_Misc.m_MotionVolumesBox.c_str())); m_Controls->m_MaxTranslationBoxX->setValue(m_Parameters.m_SignalGen.m_Translation[0]); m_Controls->m_MaxTranslationBoxY->setValue(m_Parameters.m_SignalGen.m_Translation[1]); m_Controls->m_MaxTranslationBoxZ->setValue(m_Parameters.m_SignalGen.m_Translation[2]); m_Controls->m_MaxRotationBoxX->setValue(m_Parameters.m_SignalGen.m_Rotation[0]); m_Controls->m_MaxRotationBoxY->setValue(m_Parameters.m_SignalGen.m_Rotation[1]); m_Controls->m_MaxRotationBoxZ->setValue(m_Parameters.m_SignalGen.m_Rotation[2]); m_Controls->m_Compartment1Box->setCurrentIndex(0); m_Controls->m_Compartment2Box->setCurrentIndex(0); m_Controls->m_Compartment3Box->setCurrentIndex(0); m_Controls->m_Compartment4Box->setCurrentIndex(0); for (unsigned int i=0; i* signalModel = nullptr; if (iGetVolumeFractionImage().IsNotNull() ) { compVolNode = mitk::DataNode::New(); mitk::Image::Pointer image = mitk::Image::New(); image->InitializeByItk(signalModel->GetVolumeFractionImage().GetPointer()); image->SetVolume(signalModel->GetVolumeFractionImage()->GetBufferPointer()); compVolNode->SetData( image ); compVolNode->SetName("Compartment volume "+QString::number(signalModel->m_CompartmentId).toStdString()); GetDataStorage()->Add(compVolNode); } switch (signalModel->m_CompartmentId) { case 1: { if (compVolNode.IsNotNull()) m_Controls->m_Comp1VolumeFraction->SetSelectedNode(compVolNode); if (dynamic_cast*>(signalModel)) { mitk::StickModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_StickWidget1->SetT2(model->GetT2()); m_Controls->m_StickWidget1->SetT1(model->GetT1()); m_Controls->m_StickWidget1->SetD(model->GetDiffusivity()); m_Controls->m_Compartment1Box->setCurrentIndex(0); break; } else if (dynamic_cast*>(signalModel)) { mitk::TensorModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_TensorWidget1->SetT2(model->GetT2()); m_Controls->m_TensorWidget1->SetT1(model->GetT1()); m_Controls->m_TensorWidget1->SetD1(model->GetDiffusivity1()); m_Controls->m_TensorWidget1->SetD2(model->GetDiffusivity2()); m_Controls->m_TensorWidget1->SetD3(model->GetDiffusivity3()); m_Controls->m_Compartment1Box->setCurrentIndex(2); break; } else if (dynamic_cast*>(signalModel)) { mitk::RawShModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_PrototypeWidget1->SetNumberOfSamples(model->GetMaxNumKernels()); m_Controls->m_PrototypeWidget1->SetMinFa(model->GetFaRange().first); m_Controls->m_PrototypeWidget1->SetMaxFa(model->GetFaRange().second); m_Controls->m_PrototypeWidget1->SetMinAdc(model->GetAdcRange().first); m_Controls->m_PrototypeWidget1->SetMaxAdc(model->GetAdcRange().second); m_Controls->m_Compartment1Box->setCurrentIndex(3); break; } break; } case 2: { if (compVolNode.IsNotNull()) m_Controls->m_Comp2VolumeFraction->SetSelectedNode(compVolNode); if (dynamic_cast*>(signalModel)) { mitk::StickModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_StickWidget2->SetT2(model->GetT2()); m_Controls->m_StickWidget2->SetT1(model->GetT1()); m_Controls->m_StickWidget2->SetD(model->GetDiffusivity()); m_Controls->m_Compartment2Box->setCurrentIndex(1); break; } else if (dynamic_cast*>(signalModel)) { mitk::TensorModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_TensorWidget2->SetT2(model->GetT2()); m_Controls->m_TensorWidget2->SetT1(model->GetT1()); m_Controls->m_TensorWidget2->SetD1(model->GetDiffusivity1()); m_Controls->m_TensorWidget2->SetD2(model->GetDiffusivity2()); m_Controls->m_TensorWidget2->SetD3(model->GetDiffusivity3()); m_Controls->m_Compartment2Box->setCurrentIndex(3); break; } break; } case 3: { if (compVolNode.IsNotNull()) m_Controls->m_Comp3VolumeFraction->SetSelectedNode(compVolNode); if (dynamic_cast*>(signalModel)) { mitk::BallModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_BallWidget1->SetT2(model->GetT2()); m_Controls->m_BallWidget1->SetT1(model->GetT1()); m_Controls->m_BallWidget1->SetD(model->GetDiffusivity()); m_Controls->m_Compartment3Box->setCurrentIndex(0); break; } else if (dynamic_cast*>(signalModel)) { mitk::AstroStickModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_AstrosticksWidget1->SetT2(model->GetT2()); m_Controls->m_AstrosticksWidget1->SetT1(model->GetT1()); m_Controls->m_AstrosticksWidget1->SetD(model->GetDiffusivity()); m_Controls->m_AstrosticksWidget1->SetRandomizeSticks(model->GetRandomizeSticks()); m_Controls->m_Compartment3Box->setCurrentIndex(1); break; } else if (dynamic_cast*>(signalModel)) { mitk::DotModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_DotWidget1->SetT2(model->GetT2()); m_Controls->m_DotWidget1->SetT1(model->GetT1()); m_Controls->m_Compartment3Box->setCurrentIndex(2); break; } else if (dynamic_cast*>(signalModel)) { mitk::RawShModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_PrototypeWidget3->SetNumberOfSamples(model->GetMaxNumKernels()); m_Controls->m_PrototypeWidget3->SetMinFa(model->GetFaRange().first); m_Controls->m_PrototypeWidget3->SetMaxFa(model->GetFaRange().second); m_Controls->m_PrototypeWidget3->SetMinAdc(model->GetAdcRange().first); m_Controls->m_PrototypeWidget3->SetMaxAdc(model->GetAdcRange().second); m_Controls->m_Compartment3Box->setCurrentIndex(3); break; } break; } case 4: { if (compVolNode.IsNotNull()) m_Controls->m_Comp4VolumeFraction->SetSelectedNode(compVolNode); if (dynamic_cast*>(signalModel)) { mitk::BallModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_BallWidget2->SetT2(model->GetT2()); m_Controls->m_BallWidget2->SetT1(model->GetT1()); m_Controls->m_BallWidget2->SetD(model->GetDiffusivity()); m_Controls->m_Compartment4Box->setCurrentIndex(1); break; } else if (dynamic_cast*>(signalModel)) { mitk::AstroStickModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_AstrosticksWidget2->SetT2(model->GetT2()); m_Controls->m_AstrosticksWidget2->SetT1(model->GetT1()); m_Controls->m_AstrosticksWidget2->SetD(model->GetDiffusivity()); m_Controls->m_AstrosticksWidget2->SetRandomizeSticks(model->GetRandomizeSticks()); m_Controls->m_Compartment4Box->setCurrentIndex(2); break; } else if (dynamic_cast*>(signalModel)) { mitk::DotModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_DotWidget2->SetT2(model->GetT2()); m_Controls->m_DotWidget2->SetT1(model->GetT1()); m_Controls->m_Compartment4Box->setCurrentIndex(3); break; } else if (dynamic_cast*>(signalModel)) { mitk::RawShModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_PrototypeWidget4->SetNumberOfSamples(model->GetMaxNumKernels()); m_Controls->m_PrototypeWidget4->SetMinFa(model->GetFaRange().first); m_Controls->m_PrototypeWidget4->SetMaxFa(model->GetFaRange().second); m_Controls->m_PrototypeWidget4->SetMinAdc(model->GetAdcRange().first); m_Controls->m_PrototypeWidget4->SetMaxAdc(model->GetAdcRange().second); m_Controls->m_Compartment4Box->setCurrentIndex(4); break; } break; } } } if ( m_Parameters.m_SignalGen.m_MaskImage ) { mitk::Image::Pointer image = mitk::Image::New(); image->InitializeByItk(m_Parameters.m_SignalGen.m_MaskImage.GetPointer()); image->SetVolume(m_Parameters.m_SignalGen.m_MaskImage->GetBufferPointer()); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); node->SetName("Tissue mask"); GetDataStorage()->Add(node); m_Controls->m_MaskComboBox->SetSelectedNode(node); } if ( m_Parameters.m_SignalGen.m_FrequencyMap ) { mitk::Image::Pointer image = mitk::Image::New(); image->InitializeByItk(m_Parameters.m_SignalGen.m_FrequencyMap.GetPointer()); image->SetVolume(m_Parameters.m_SignalGen.m_FrequencyMap->GetBufferPointer()); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); node->SetName("Frequency map"); GetDataStorage()->Add(node); m_Controls->m_FrequencyMapBox->SetSelectedNode(node); } } void QmitkFiberfoxView::ShowAdvancedOptions(int state) { if (state) { m_Controls->m_AdvancedSignalOptionsFrame->setVisible(true); m_Controls->m_AdvancedOptionsBox_2->setChecked(true); } else { m_Controls->m_AdvancedSignalOptionsFrame->setVisible(false); m_Controls->m_AdvancedOptionsBox_2->setChecked(false); } } void QmitkFiberfoxView::Comp1ModelFrameVisibility(int index) { m_Controls->m_StickWidget1->setVisible(false); m_Controls->m_ZeppelinWidget1->setVisible(false); m_Controls->m_TensorWidget1->setVisible(false); m_Controls->m_PrototypeWidget1->setVisible(false); switch (index) { case 0: m_Controls->m_StickWidget1->setVisible(true); break; case 1: m_Controls->m_ZeppelinWidget1->setVisible(true); break; case 2: m_Controls->m_TensorWidget1->setVisible(true); break; case 3: m_Controls->m_PrototypeWidget1->setVisible(true); break; } } void QmitkFiberfoxView::Comp2ModelFrameVisibility(int index) { m_Controls->m_StickWidget2->setVisible(false); m_Controls->m_ZeppelinWidget2->setVisible(false); m_Controls->m_TensorWidget2->setVisible(false); m_Controls->m_Comp2FractionFrame->setVisible(false); switch (index) { case 0: break; case 1: m_Controls->m_StickWidget2->setVisible(true); m_Controls->m_Comp2FractionFrame->setVisible(true); break; case 2: m_Controls->m_ZeppelinWidget2->setVisible(true); m_Controls->m_Comp2FractionFrame->setVisible(true); break; case 3: m_Controls->m_TensorWidget2->setVisible(true); m_Controls->m_Comp2FractionFrame->setVisible(true); break; } } void QmitkFiberfoxView::Comp3ModelFrameVisibility(int index) { m_Controls->m_BallWidget1->setVisible(false); m_Controls->m_AstrosticksWidget1->setVisible(false); m_Controls->m_DotWidget1->setVisible(false); m_Controls->m_PrototypeWidget3->setVisible(false); switch (index) { case 0: m_Controls->m_BallWidget1->setVisible(true); break; case 1: m_Controls->m_AstrosticksWidget1->setVisible(true); break; case 2: m_Controls->m_DotWidget1->setVisible(true); break; case 3: m_Controls->m_PrototypeWidget3->setVisible(true); break; } } void QmitkFiberfoxView::Comp4ModelFrameVisibility(int index) { m_Controls->m_BallWidget2->setVisible(false); m_Controls->m_AstrosticksWidget2->setVisible(false); m_Controls->m_DotWidget2->setVisible(false); m_Controls->m_PrototypeWidget4->setVisible(false); m_Controls->m_Comp4FractionFrame->setVisible(false); switch (index) { case 0: break; case 1: m_Controls->m_BallWidget2->setVisible(true); m_Controls->m_Comp4FractionFrame->setVisible(true); break; case 2: m_Controls->m_AstrosticksWidget2->setVisible(true); m_Controls->m_Comp4FractionFrame->setVisible(true); break; case 3: m_Controls->m_DotWidget2->setVisible(true); m_Controls->m_Comp4FractionFrame->setVisible(true); break; case 4: m_Controls->m_PrototypeWidget4->setVisible(true); m_Controls->m_Comp4FractionFrame->setVisible(true); break; } } void QmitkFiberfoxView::OnAddMotion(int value) { if (value>0) m_Controls->m_MotionArtifactFrame->setVisible(true); else m_Controls->m_MotionArtifactFrame->setVisible(false); } void QmitkFiberfoxView::OnAddDrift(int value) { if (value>0) m_Controls->m_DriftFrame->setVisible(true); else m_Controls->m_DriftFrame->setVisible(false); } void QmitkFiberfoxView::OnAddAliasing(int value) { if (value>0) m_Controls->m_AliasingFrame->setVisible(true); else m_Controls->m_AliasingFrame->setVisible(false); } void QmitkFiberfoxView::OnAddSpikes(int value) { if (value>0) m_Controls->m_SpikeFrame->setVisible(true); else m_Controls->m_SpikeFrame->setVisible(false); } void QmitkFiberfoxView::OnAddEddy(int value) { if (value>0) m_Controls->m_EddyFrame->setVisible(true); else m_Controls->m_EddyFrame->setVisible(false); } void QmitkFiberfoxView::OnAddDistortions(int value) { if (value>0) m_Controls->m_DistortionsFrame->setVisible(true); else m_Controls->m_DistortionsFrame->setVisible(false); } void QmitkFiberfoxView::OnAddGhosts(int value) { if (value>0) m_Controls->m_GhostFrame->setVisible(true); else m_Controls->m_GhostFrame->setVisible(false); } void QmitkFiberfoxView::OnAddNoise(int value) { if (value>0) m_Controls->m_NoiseFrame->setVisible(true); else m_Controls->m_NoiseFrame->setVisible(false); } +void QmitkFiberfoxView::OnAddRinging(int value) +{ + if (value>0) + m_Controls->m_ZeroRinging->setVisible(true); + else + m_Controls->m_ZeroRinging->setVisible(false); +} + QmitkFiberfoxView::GradientListType QmitkFiberfoxView::GenerateHalfShell(int NPoints) { NPoints *= 2; GradientListType pointshell; int numB0 = NPoints/20; if (numB0==0) numB0=1; GradientType g; g.Fill(0.0); for (int i=0; i theta; theta.set_size(NPoints); vnl_vector phi; phi.set_size(NPoints); double C = sqrt(4*itk::Math::pi); phi(0) = 0.0; phi(NPoints-1) = 0.0; for(int i=0; i0 && i std::vector > QmitkFiberfoxView::MakeGradientList() { std::vector > retval; vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); // Add 0 vector for B0 int numB0 = ndirs/10; if (numB0==0) numB0=1; itk::Vector v; v.Fill(0.0); for (int i=0; i v; v[0] = U->get(0,i); v[1] = U->get(1,i); v[2] = U->get(2,i); retval.push_back(v); } return retval; } void QmitkFiberfoxView::GenerateImage() { if (m_Controls->m_FiberBundleComboBox->GetSelectedNode().IsNull() && !mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( m_Controls->m_TemplateComboBox->GetSelectedNode())) { mitk::Image::Pointer image = mitk::ImageGenerator::GenerateGradientImage( m_Controls->m_SizeX->value(), m_Controls->m_SizeY->value(), m_Controls->m_SizeZ->value(), m_Controls->m_SpacingX->value(), m_Controls->m_SpacingY->value(), m_Controls->m_SpacingZ->value()); mitk::Point3D origin; origin[0] = m_Controls->m_SpacingX->value()/2; origin[1] = m_Controls->m_SpacingY->value()/2; origin[2] = m_Controls->m_SpacingZ->value()/2; image->SetOrigin(origin); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); node->SetName("Dummy"); unsigned int window = m_Controls->m_SizeX->value()*m_Controls->m_SizeY->value()*m_Controls->m_SizeZ->value(); unsigned int level = window/2; mitk::LevelWindow lw; lw.SetLevelWindow(level, window); node->SetProperty( "levelwindow", mitk::LevelWindowProperty::New( lw ) ); GetDataStorage()->Add(node); m_SelectedImageNode = node; mitk::BaseData::Pointer basedata = node->GetData(); if (basedata.IsNotNull()) { mitk::RenderingManager::GetInstance()->InitializeViews( basedata->GetTimeGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true ); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } UpdateGui(); QMessageBox::information(nullptr, "Template image generated", "You have selected no fiber bundle or diffusion-weighted image, which can be used to simulate a new diffusion-weighted image. A template image with the specified geometry has been generated that can be used to draw artificial fibers (see view 'Fiber Generator')."); } else if (m_Controls->m_FiberBundleComboBox->GetSelectedNode().IsNotNull()) SimulateImageFromFibers(m_Controls->m_FiberBundleComboBox->GetSelectedNode()); else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( m_Controls->m_TemplateComboBox->GetSelectedNode()) ) SimulateForExistingDwi(m_Controls->m_TemplateComboBox->GetSelectedNode()); else QMessageBox::information(nullptr, "No image generated", "You have selected no fiber bundle or diffusion-weighted image, which can be used to simulate a new diffusion-weighted image."); } void QmitkFiberfoxView::SetFocus() { } void QmitkFiberfoxView::SimulateForExistingDwi(mitk::DataNode* imageNode) { bool isDiffusionImage( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(imageNode->GetData())) ); if ( !isDiffusionImage ) { return; } UpdateParametersFromGui(); mitk::Image::Pointer diffImg = dynamic_cast(imageNode->GetData()); ItkDwiType::Pointer itkVectorImagePointer = ItkDwiType::New(); mitk::CastToItkImage(diffImg, itkVectorImagePointer); m_TractsToDwiFilter = itk::TractsToDWIImageFilter< short >::New(); m_Parameters.m_Misc.m_ParentNode = imageNode; m_Parameters.m_SignalGen.m_SignalScale = 1; m_Parameters.m_Misc.m_ResultNode->SetName(m_Parameters.m_Misc.m_ParentNode->GetName() +"_D"+QString::number(m_Parameters.m_SignalGen.m_ImageRegion.GetSize(0)).toStdString() +"-"+QString::number(m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1)).toStdString() +"-"+QString::number(m_Parameters.m_SignalGen.m_ImageRegion.GetSize(2)).toStdString() +"_S"+QString::number(m_Parameters.m_SignalGen.m_ImageSpacing[0]).toStdString() +"-"+QString::number(m_Parameters.m_SignalGen.m_ImageSpacing[1]).toStdString() +"-"+QString::number(m_Parameters.m_SignalGen.m_ImageSpacing[2]).toStdString() +"_b"+QString::number(m_Parameters.m_SignalGen.GetBvalue()).toStdString() +"_"+m_Parameters.m_Misc.m_SignalModelString +m_Parameters.m_Misc.m_ArtifactModelString); m_Parameters.ApplyDirectionMatrix(); m_TractsToDwiFilter->SetParameters(m_Parameters); m_TractsToDwiFilter->SetInputImage(itkVectorImagePointer); m_Thread.start(QThread::LowestPriority); } void QmitkFiberfoxView::SimulateImageFromFibers(mitk::DataNode* fiberNode) { mitk::FiberBundle::Pointer fiberBundle = dynamic_cast(fiberNode->GetData()); if (fiberBundle->GetNumFibers()<=0) { return; } UpdateParametersFromGui(); m_TractsToDwiFilter = itk::TractsToDWIImageFilter< short >::New(); m_Parameters.m_Misc.m_ParentNode = fiberNode; m_Parameters.m_Misc.m_ResultNode->SetName(m_Parameters.m_Misc.m_ParentNode->GetName() +"_D"+QString::number(m_Parameters.m_SignalGen.m_ImageRegion.GetSize(0)).toStdString() +"-"+QString::number(m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1)).toStdString() +"-"+QString::number(m_Parameters.m_SignalGen.m_ImageRegion.GetSize(2)).toStdString() +"_S"+QString::number(m_Parameters.m_SignalGen.m_ImageSpacing[0]).toStdString() +"-"+QString::number(m_Parameters.m_SignalGen.m_ImageSpacing[1]).toStdString() +"-"+QString::number(m_Parameters.m_SignalGen.m_ImageSpacing[2]).toStdString() +"_b"+QString::number(m_Parameters.m_SignalGen.GetBvalue()).toStdString() +"_"+m_Parameters.m_Misc.m_SignalModelString +m_Parameters.m_Misc.m_ArtifactModelString); if ( m_Controls->m_TemplateComboBox->GetSelectedNode().IsNotNull() && mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast (m_Controls->m_TemplateComboBox->GetSelectedNode()->GetData()) ) ) { bool first = true; bool ok = true; mitk::Image::Pointer diffImg = dynamic_cast(m_Controls->m_TemplateComboBox->GetSelectedNode()->GetData()); itk::Image< itk::DiffusionTensor3D< double >, 3 >::Pointer tensorImage = nullptr; const int shOrder = 2; typedef itk::AnalyticalDiffusionQballReconstructionImageFilter QballFilterType; QballFilterType::CoefficientImageType::Pointer itkFeatureImage = nullptr; ItkDoubleImgType::Pointer adcImage = nullptr; for (unsigned int i=0; i* model = nullptr; if (i* >(m_Parameters.m_FiberModelList.at(i)); else model = dynamic_cast< mitk::RawShModel<>* >(m_Parameters.m_NonFiberModelList.at(i-m_Parameters.m_FiberModelList.size())); if (model!=0 && model->GetNumberOfKernels()<=0) { if (first==true) { ItkDwiType::Pointer itkVectorImagePointer = ItkDwiType::New(); mitk::CastToItkImage(diffImg, itkVectorImagePointer); typedef itk::DiffusionTensor3DReconstructionImageFilter< short, short, double > TensorReconstructionImageFilterType; TensorReconstructionImageFilterType::Pointer filter = TensorReconstructionImageFilterType::New(); filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(diffImg)); filter->SetGradientImage(mitk::DiffusionPropertyHelper::GetGradientContainer(diffImg), itkVectorImagePointer ); filter->Update(); tensorImage = filter->GetOutput(); QballFilterType::Pointer qballfilter = QballFilterType::New(); qballfilter->SetGradientImage(mitk::DiffusionPropertyHelper::GetGradientContainer(diffImg), itkVectorImagePointer ); qballfilter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(diffImg)); qballfilter->SetLambda(0.006); qballfilter->SetNormalizationMethod(QballFilterType::QBAR_RAW_SIGNAL); qballfilter->Update(); itkFeatureImage = qballfilter->GetCoefficientImage(); itk::AdcImageFilter< short, double >::Pointer adcFilter = itk::AdcImageFilter< short, double >::New(); adcFilter->SetInput( itkVectorImagePointer ); adcFilter->SetGradientDirections(mitk::DiffusionPropertyHelper::GetGradientContainer(diffImg)); adcFilter->SetB_value(mitk::DiffusionPropertyHelper::GetReferenceBValue(diffImg)); adcFilter->Update(); adcImage = adcFilter->GetOutput(); } ok = model->SampleKernels(diffImg, m_Parameters.m_SignalGen.m_MaskImage, tensorImage, itkFeatureImage, adcImage); if (!ok) break; } } if (!ok) { QMessageBox::information( nullptr, "Simulation cancelled", "No valid prototype signals could be sampled."); return; } } else if ( m_Controls->m_Compartment1Box->currentIndex()==3 || m_Controls->m_Compartment3Box->currentIndex()==3 || m_Controls->m_Compartment4Box->currentIndex()==4 ) { QMessageBox::information( nullptr, "Simulation cancelled", "Prototype signal but no diffusion-weighted image selected to sample signal from."); return; } m_Parameters.ApplyDirectionMatrix(); m_TractsToDwiFilter->SetParameters(m_Parameters); m_TractsToDwiFilter->SetFiberBundle(fiberBundle); m_Thread.start(QThread::LowestPriority); } void QmitkFiberfoxView::SetBvalsEdit() { // SELECT FOLDER DIALOG std::string filename; filename = QFileDialog::getOpenFileName(nullptr, "Select bvals file", QString(filename.c_str())).toStdString(); if (filename.empty()) m_Controls->m_LoadBvalsEdit->setText("-"); else m_Controls->m_LoadBvalsEdit->setText(QString(filename.c_str())); } void QmitkFiberfoxView::SetBvecsEdit() { // SELECT FOLDER DIALOG std::string filename; filename = QFileDialog::getOpenFileName(nullptr, "Select bvecs file", QString(filename.c_str())).toStdString(); if (filename.empty()) m_Controls->m_LoadBvecsEdit->setText("-"); else m_Controls->m_LoadBvecsEdit->setText(QString(filename.c_str())); } void QmitkFiberfoxView::SetOutputPath() { // SELECT FOLDER DIALOG std::string outputPath; outputPath = QFileDialog::getExistingDirectory(nullptr, "Save images to...", QString(outputPath.c_str())).toStdString(); if (outputPath.empty()) m_Controls->m_SavePathEdit->setText("-"); else { outputPath += "/"; m_Controls->m_SavePathEdit->setText(QString(outputPath.c_str())); } } void QmitkFiberfoxView::UpdateGui() { m_Controls->m_GeometryFrame->setEnabled(true); m_Controls->m_GeometryMessage->setVisible(false); m_Controls->m_DiffusionPropsMessage->setVisible(false); m_Controls->m_LoadGradientsFrame->setVisible(false); m_Controls->m_GenerateGradientsFrame->setVisible(false); if (m_Controls->m_UseBvalsBvecsBox->isChecked()) m_Controls->m_LoadGradientsFrame->setVisible(true); else m_Controls->m_GenerateGradientsFrame->setVisible(true); // Signal generation gui if (m_Controls->m_MaskComboBox->GetSelectedNode().IsNotNull() || m_Controls->m_TemplateComboBox->GetSelectedNode().IsNotNull()) { m_Controls->m_GeometryMessage->setVisible(true); m_Controls->m_GeometryFrame->setEnabled(false); } if ( m_Controls->m_TemplateComboBox->GetSelectedNode().IsNotNull() && mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_Controls->m_TemplateComboBox->GetSelectedNode()->GetData()) ) ) { m_Controls->m_DiffusionPropsMessage->setVisible(true); m_Controls->m_GeometryMessage->setVisible(true); m_Controls->m_GeometryFrame->setEnabled(false); m_Controls->m_LoadGradientsFrame->setVisible(false); m_Controls->m_GenerateGradientsFrame->setVisible(false); } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberfox/src/internal/QmitkFiberfoxView.h b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberfox/src/internal/QmitkFiberfoxView.h index 3d9ea22d68..092bd66e4f 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberfox/src/internal/QmitkFiberfoxView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberfox/src/internal/QmitkFiberfoxView.h @@ -1,171 +1,172 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include "ui_QmitkFiberfoxViewControls.h" #include #include #include #ifndef Q_MOC_RUN #include #include #include #include #include #include #include #include #include #include #include #include #include #include #endif #include #include #include #include /*! \brief View for fiber based diffusion software phantoms (Fiberfox). See "Fiberfox: Facilitating the creation of realistic white matter software phantoms" (DOI: 10.1002/mrm.25045) for details. */ // Forward Qt class declarations class QmitkFiberfoxView; class QmitkFiberfoxWorker : public QObject { Q_OBJECT public: QmitkFiberfoxWorker(QmitkFiberfoxView* view); public slots: void run(); private: QmitkFiberfoxView* m_View; }; class QmitkFiberfoxView : public QmitkAbstractView { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: static const std::string VIEW_ID; QmitkFiberfoxView(); virtual ~QmitkFiberfoxView(); virtual void CreateQtPartControl(QWidget *parent) override; void SetFocus() override; typedef mitk::DiffusionPropertyHelper::GradientDirectionType GradientDirectionType; typedef mitk::DiffusionPropertyHelper::GradientDirectionsContainerType GradientDirectionContainerType; typedef itk::Vector GradientType; typedef std::vector GradientListType; typedef itk::VectorImage< short, 3 > ItkDwiType; typedef itk::Image ItkDoubleImgType; typedef itk::Image ItkFloatImgType; typedef itk::Image ItkUcharImgType; template std::vector > MakeGradientList(); protected slots: void SetOutputPath(); ///< path where image is automatically saved to after the simulation is finished void LoadParameters(); ///< load fiberfox parameters void SaveParameters(); ///< save fiberfox parameters void SetBvalsEdit(); void SetBvecsEdit(); void BeforeThread(); void AfterThread(); void KillThread(); ///< abort simulation void UpdateSimulationStatus(); ///< print simulation progress and satus messages void GenerateImage(); ///< start image simulation void Comp1ModelFrameVisibility(int index); ///< only show parameters of selected signal model for compartment 1 void Comp2ModelFrameVisibility(int index); ///< only show parameters of selected signal model for compartment 2 void Comp3ModelFrameVisibility(int index); ///< only show parameters of selected signal model for compartment 3 void Comp4ModelFrameVisibility(int index); ///< only show parameters of selected signal model for compartment 4 void ShowAdvancedOptions(int state); /** update GUI elements */ void OnAddNoise(int value); void OnAddGhosts(int value); void OnAddDistortions(int value); void OnAddEddy(int value); void OnAddSpikes(int value); void OnAddAliasing(int value); void OnAddMotion(int value); void OnAddDrift(int value); + void OnAddRinging(int value); void OnMaskSelected(int value); void OnFibSelected(int value); void OnTemplateSelected(int value); void OnBvalsBvecsCheck(int value); protected: GradientListType GenerateHalfShell(int NPoints); ///< generate vectors distributed over the halfsphere Ui::QmitkFiberfoxViewControls* m_Controls; void SimulateForExistingDwi(mitk::DataNode* imageNode); ///< add artifacts to existing diffusion weighted image void SimulateImageFromFibers(mitk::DataNode* fiberNode); ///< simulate new diffusion weighted image void UpdateParametersFromGui(); ///< update fiberfox paramater object void UpdateGui(); ///< enable/disbale buttons etc. according to current datamanager selection void PlanarFigureSelected( itk::Object* object, const itk::EventObject& ); void EnableCrosshairNavigation(); ///< enable crosshair navigation if planar figure interaction ends void DisableCrosshairNavigation(); ///< disable crosshair navigation if planar figure interaction starts void SaveParameters(QString filename); mitk::DataNode::Pointer m_SelectedImageNode; QString m_ParameterFile; ///< parameter file name // GUI thread QmitkFiberfoxWorker m_Worker; ///< runs filter QThread m_Thread; ///< worker thread bool m_ThreadIsRunning; QTimer* m_SimulationTimer; QTime m_SimulationTime; QString m_SimulationStatusText; /** Image filters that do all the simulations. */ itk::TractsToDWIImageFilter< short >::Pointer m_TractsToDwiFilter; friend class QmitkFiberfoxWorker; FiberfoxParameters m_Parameters; }; 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 c9f768a18c..0548f27733 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,2742 +1,2755 @@ QmitkFiberfoxViewControls 0 0 463 - 2651 + 2713 Form Save Parameters :/QmitkDiffusionImaging/general_icons/download.ico:/QmitkDiffusionImaging/general_icons/download.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> Start Simulation :/QmitkDiffusionImaging/general_icons/right.ico:/QmitkDiffusionImaging/general_icons/right.ico QGroupBox { background-color: transparent; } Input Data 6 6 6 6 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 - ... <html><head/><body><p>Select a binary image to define the area of signal generation. Outside of the mask image only noise will be actively generated.</p></body></html> QComboBox::AdjustToMinimumContentsLength Fiber Bundle: false Save path: false Tissue Mask: false <html><head/><body><p>Select a fiber bundle to generate the white matter signal from. You can either use the fiber definition tab to manually define an input fiber bundle or you can also use any existing bundle, e.g. yielded by a tractography algorithm.</p></body></html> QComboBox::AdjustToMinimumContentsLength Template Image: false <html><head/><body><p>The parameters for the simulation (e.g. spacing, size, diffuison-weighted gradients, b-value) are adopted from this image.</p></body></html> QComboBox::AdjustToMinimumContentsLength 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 6 6 6 6 Select signal model for intra-axonal compartment. Stick Model Zeppelin Model Tensor Model Prototype Signal QFrame::NoFrame QFrame::Raised 0 0 0 0 Volume Fraction: Optional! If no volume fraction map for this compartment is set, the corresponding volume fractions are calculated from the input fibers. Load Parameters :/QmitkDiffusionImaging/general_icons/upload.ico:/QmitkDiffusionImaging/general_icons/upload.ico QGroupBox { background-color: transparent; } Extra-axonal Compartments 6 6 6 6 QFrame::NoFrame QFrame::Raised 0 0 0 0 Volume Fraction: Select signal model for extra-axonal compartment. Ball Model Astrosticks Model Dot Model Prototype Signal Qt::Horizontal QFrame::NoFrame QFrame::Raised 0 0 0 0 Volume Fraction: Optional! If no volume fraction map for this compartment is set, the corresponding volume fractions are calculated from the input fibers. Select signal model for extra-axonal compartment. -- Ball Model Astrosticks Model Dot Model Prototype Signal QGroupBox { background-color: transparent; } Noise and other Artifacts 6 6 6 6 - - - - Add Distortions - - - false - - - - - - - Add ringing artifacts occuring at strong edges in the image. - - - Add Gibbs Ringing - - - false - - - - - - - Add Motion Artifacts - - - false - - - - - - - Add Aliasing - - - false - - - - - - - Qt::Horizontal - - - - - - - Add Signal Drift - - - false - - - - - - - QFrame::NoFrame - - - QFrame::Raised - - - - 0 - - - 0 - - - 0 - - - 0 - - - - - Num. Spikes: - - - - - - - The number of randomly occurring signal spikes. - - - 1 - - - - - - - Spike amplitude relative to the largest signal amplitude of the corresponding k-space slice. - - - 0.100000000000000 - - - 0.100000000000000 - - - - - - - Scale: - - - - - - - - - - Qt::Horizontal - - - - - - - Qt::Horizontal - - - - - - - Qt::Horizontal - - - - - - - true - - - QFrame::NoFrame - - - QFrame::Raised - - - - 6 - - - 0 - - - 0 - - - 0 - - - 0 - - - - - - - - - - - - - - Shrink FOV (%): - - - false - - - - - - - Shrink FOV by this percentage. - - - 1 - - - 0.000000000000000 - - - 90.000000000000000 - - - 0.100000000000000 - - - 25.000000000000000 - - - - - - - - - - true - - - QFrame::NoFrame - - - QFrame::Raised - - - - 6 - - - 0 - - - 0 - - - 0 - - - 0 - - - - - - - - - - - - - - Frequency Map: - - - false - - - - - - - Select image specifying the frequency inhomogeneities (in Hz). - - - - - - - - - - Qt::Horizontal - - - - - - - Qt::Horizontal - - - - - - - QFrame::NoFrame - - - QFrame::Raised - - - - 0 - - - 0 - - - 0 - - - 0 - - - - - Variance: - - - - - - - Variance of selected noise distribution. - - - 10 - - - 0.000000000000000 - - - 999999999.000000000000000 - - - 0.001000000000000 - - - 50.000000000000000 - - - - - - - Distribution: - - - - - - - Noise distribution - - - - Complex Gaussian - - - - - Rician - - - - - - - - - - - - - - Add Eddy Current Effects - - - false - - - - - - - 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 - - - - + + + + Qt::Horizontal + true QFrame::NoFrame QFrame::Raised 0 6 0 0 6 Toggle between random movement and linear movement. Randomize motion true QGroupBox { background-color: transparent; } 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 QGroupBox { background-color: transparent; } 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 QFrame::NoFrame QFrame::Raised 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 - - - - Add Noise + + + + QFrame::NoFrame - - false + + QFrame::Raised + + + 0 + + + 0 + + + 0 + + + 0 + + + + + Num. Spikes: + + + + + + + The number of randomly occurring signal spikes. + + + 1 + + + + + + + Spike amplitude relative to the largest signal amplitude of the corresponding k-space slice. + + + 0.100000000000000 + + + 0.100000000000000 + + + + + + + Scale: + + + + - - - - Add Spikes + + + + true - - false + + QFrame::NoFrame + + + QFrame::Raised + + + + 6 + + + 0 + + + 0 + + + 0 + + + 0 + + + + + + + + + + + + + + Shrink FOV (%): + + + false + + + + + + + Shrink FOV by this percentage. + + + 1 + + + 0.000000000000000 + + + 90.000000000000000 + + + 0.100000000000000 + + + 25.000000000000000 + + + + + + + + + + Qt::Horizontal + + + + + + + true + + + QFrame::NoFrame + + + QFrame::Raised + + + + 6 + + + 0 + + + 0 + + + 0 + + + 0 + + + + + + + + + + + + + + Signal Reduction (%): + + + false + + + + + + + Global signal in last simulated volume is specified percentage lower than in the first volume. + + + 1 + + + 100.000000000000000 + + + 1.000000000000000 + + + 6.000000000000000 + + + + + + + + + + true + + QFrame::NoFrame + + + QFrame::Raised + + + + 6 + + + 0 + + + 0 + + + 0 + + + 0 + + + + + + + + + + + + + + Frequency Map: + + + false + + + + + + + Select image specifying the frequency inhomogeneities (in Hz). + + + + true QFrame::NoFrame QFrame::Raised QFormLayout::AllNonFixedFieldsGrow 6 0 0 0 0 Gradient: false Eddy current induced magnetic field gradient (in mT/m). 4 1000.000000000000000 0.001000000000000 0.010000000000000 - - + + + + Qt::Horizontal + + + + + + + + - Add N/2 Ghosts + Add Eddy Current Effects false - - + + + + Add Distortions + + + false + + + + + + + Add Spikes + + + false + + + + + + + Add Signal Drift + + + false + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + 0 + + + 0 + + + 0 + + + + + Variance: + + + + + + + Variance of selected noise distribution. + + + 10 + + + 0.000000000000000 + + + 999999999.000000000000000 + + + 0.001000000000000 + + + 50.000000000000000 + + + + + + + Distribution: + + + + + + + Noise distribution + + + + Complex Gaussian + + + + + Rician + + + + + + + + + + + Qt::Horizontal + + + + + + + Qt::Horizontal + + + + + true QFrame::NoFrame QFrame::Raised - + 6 0 0 0 0 - + - Signal Reduction (%): + K-Space Line Offset: false - + - Global signal in last simulated volume is specified percentage lower than in the first volume. + A larger offset increases the inensity of the ghost image. - 1 + 3 - 100.000000000000000 + 1.000000000000000 - 1.000000000000000 + 0.010000000000000 - 6.000000000000000 + 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! QFrame::NoFrame QFrame::Raised 0 0 0 0 6 <html><head/><body><p>Number of Channels:</p></body></html> false TR in milliseconds 1 999999999 1 4000 Output one image per compartment containing the corresponding volume fractions per voxel. Reverse Phase Encoding Direction false Signal Scale: Dwell time (time to read one line in k-space) in ms. 100.000000000000000 0.100000000000000 1.000000000000000 TE in milliseconds 1 999999999 1 100 Fiber Radius: Partial Fourier: false 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>Coil Sensitivity:</p></body></html> false Partial fourier factor (0.5-1) 3 0.500000000000000 1.000000000000000 0.100000000000000 1.000000000000000 Output phase image and volume fraction maps. Output Additional Images false Relaxation time due to magnetic field inhomogeneities (T2', in milliseconds). 1 10000 1 50 <html><head/><body><p>Repetition Time <span style=" font-style:italic;">TR</span>: </p></body></html> false Constant Linear Exponential Dwell Time: false <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 TE in milliseconds 1 10000 1 100 Number of coil elements used for the acquisiton. 1 128 1 1 Acquisition Type: Single Shot EPI Spin Echo Fiber radius used to calculate volume fractions (in µm). Set to 0 for automatic radius estimation. 9999.000000000000000 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 Gradient Directions: 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 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 1000 1 20 Fiber sampling factor which determines the accuracy of the calculated fiber and non-fiber volume fractions. 1 1000 1 3 Use bvals/bvecs files QFrame::NoFrame QFrame::Raised 0 0 0 0 ... ... - false Bvecs: false Bvals: false - false QGroupBox { background-color: transparent; } Inter-axonal Compartment 6 6 6 6 Select signal model for intra-axonal compartment. -- Stick Model Zeppelin Model Tensor Model QFrame::NoFrame QFrame::Raised 0 0 0 0 Volume Fraction: Optional! If no volume fraction map for this compartment is set, the corresponding volume fractions are calculated from the input fibers. 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_SizeX m_SizeY m_SizeZ m_SpacingX m_SpacingY m_SpacingZ m_NumGradientsBox m_BvalueBox m_AdvancedOptionsBox_2 m_SignalScaleBox m_TEbox m_TRbox 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_AddMotion m_RandomMotion m_MaxRotationBoxX m_MaxRotationBoxY m_MaxRotationBoxZ m_MaxTranslationBoxX m_MaxTranslationBoxY m_MaxTranslationBoxZ m_AddEddy m_EddyGradientStrength m_AddGibbsRinging m_SaveParametersButton m_LoadParametersButton