diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkDftImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkDftImageFilter.cpp index 89a28b3a0d..1c58e38e3e 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkDftImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkDftImageFilter.cpp @@ -1,97 +1,119 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkDftImageFilter_txx #define __itkDftImageFilter_txx #include #include #include #include "itkDftImageFilter.h" #include #include #include #define _USE_MATH_DEFINES #include namespace itk { template< class TPixelType > DftImageFilter< TPixelType > ::DftImageFilter() { this->SetNumberOfRequiredInputs( 1 ); } template< class TPixelType > void DftImageFilter< TPixelType > ::BeforeThreadedGenerateData() { } template< class TPixelType > void DftImageFilter< TPixelType > ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType) { typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); typedef ImageRegionConstIterator< InputImageType > InputIteratorType; typename InputImageType::Pointer inputImage = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); int szx = outputImage->GetLargestPossibleRegion().GetSize(0); int szy = outputImage->GetLargestPossibleRegion().GetSize(1); while( !oit.IsAtEnd() ) { - int kx = oit.GetIndex()[0]; - int ky = oit.GetIndex()[1]; + double kx = oit.GetIndex()[0]; + double ky = oit.GetIndex()[1]; +// kx -= (double)szx/2; +// ky -= (double)szy/2; - if( kx < szx/2 ) - kx = kx + szx/2; + if ((int)szx%2==1) + kx -= (szx-1)/2; else - kx = kx - szx/2; - - if( ky < szy/2 ) - ky = ky + szy/2; + kx -= szx/2; + if ((int)szy%2==1) + ky -= (szy-1)/2; else - ky = ky - szy/2; + ky -= szy/2; + +// if( kx < szx/2 ) +// kx = kx + szx/2; +// else +// kx = kx - szx/2; + +// if( ky < szy/2 ) +// ky = ky + szy/2; +// else +// ky = ky - szy/2; vcl_complex s(0,0); InputIteratorType it(inputImage, inputImage->GetLargestPossibleRegion() ); while( !it.IsAtEnd() ) { int x = it.GetIndex()[0]; int y = it.GetIndex()[1]; +// x -= (double)szx/2; +// y -= (double)szy/2; + + if ((int)szx%2==1) + x -= (szx-1)/2; + else + x -= szx/2; + if ((int)szy%2==1) + y -= (szy-1)/2; + else + y -= szy/2; vcl_complex f(it.Get().real(), it.Get().imag()); s += f * exp( std::complex(0, -2 * M_PI * (kx*(double)x/szx + ky*(double)y/szy) ) ); ++it; } oit.Set(s); ++oit; } } } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.cpp index 929497fbc7..6c8c0bd12b 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.cpp @@ -1,384 +1,373 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkKspaceImageFilter_txx #define __itkKspaceImageFilter_txx #include #include #include #include "itkKspaceImageFilter.h" #include #include #include #include #define _USE_MATH_DEFINES #include namespace itk { template< class TPixelType > KspaceImageFilter< TPixelType > ::KspaceImageFilter() : m_Z(0) , m_UseConstantRandSeed(false) , m_SpikesPerSlice(0) , m_IsBaseline(true) { m_DiffusionGradientDirection.Fill(0.0); m_CoilPosition.Fill(0.0); } template< class TPixelType > void KspaceImageFilter< TPixelType > ::BeforeThreadedGenerateData() { m_Spike = vcl_complex(0,0); typename OutputImageType::Pointer outputImage = OutputImageType::New(); itk::ImageRegion<2> region; region.SetSize(0, m_OutSize[0]); region.SetSize(1, m_OutSize[1]); outputImage->SetLargestPossibleRegion( region ); outputImage->SetBufferedRegion( region ); outputImage->SetRequestedRegion( region ); outputImage->Allocate(); outputImage->FillBuffer(m_Spike); m_KSpaceImage = InputImageType::New(); m_KSpaceImage->SetLargestPossibleRegion( region ); m_KSpaceImage->SetBufferedRegion( region ); m_KSpaceImage->SetRequestedRegion( region ); m_KSpaceImage->Allocate(); m_KSpaceImage->FillBuffer(0.0); m_Gamma = 42576000; // Gyromagnetic ratio in Hz/T if ( m_Parameters.m_SignalGen.m_EddyStrength>0 && m_DiffusionGradientDirection.GetNorm()>0.001) { m_DiffusionGradientDirection.Normalize(); m_DiffusionGradientDirection = m_DiffusionGradientDirection * m_Parameters.m_SignalGen.m_EddyStrength/1000 * m_Gamma; m_IsBaseline = false; } this->SetNthOutput(0, outputImage); m_Transform = m_Parameters.m_SignalGen.m_ImageDirection; for (int i=0; i<3; i++) for (int j=0; j<3; j++) m_Transform[i][j] *= m_Parameters.m_SignalGen.m_ImageSpacing[j]; 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 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; } } } template< class TPixelType > double KspaceImageFilter< TPixelType >::CoilSensitivity(DoubleVectorType& 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: { DoubleVectorType diff = pos-m_CoilPosition; double sens = diff.GetNorm()*m_CoilSensitivityFactor + 1; if (sens<0) sens = 0; return sens; } case SignalGenerationParameters::COIL_EXPONENTIAL: { DoubleVectorType diff = pos-m_CoilPosition; double dist = diff.GetNorm(); return exp(-dist*m_CoilSensitivityFactor); } default: return 1; } } template< class TPixelType > void KspaceImageFilter< TPixelType > ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType threadID) { itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer randGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); randGen->SetSeed(); if (m_UseConstantRandSeed) // always generate the same random numbers? randGen->SetSeed(threadID*100); else randGen->SetSeed(); typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); typedef ImageRegionConstIterator< InputImageType > InputIteratorType; double kxMax = outputImage->GetLargestPossibleRegion().GetSize(0); // k-space size in x-direction double kyMax = outputImage->GetLargestPossibleRegion().GetSize(1); // k-space size in y-direction double xMax = m_CompartmentImages.at(0)->GetLargestPossibleRegion().GetSize(0); // scanner coverage in x-direction double yMax = m_CompartmentImages.at(0)->GetLargestPossibleRegion().GetSize(1); // scanner coverage in y-direction double numPix = kxMax*kyMax; double dt = m_Parameters.m_SignalGen.m_tLine/kxMax; - double fromMaxEcho = -m_Parameters.m_SignalGen.m_tLine*kyMax/2-dt*kxMax/2; - double upsampling = xMax/kxMax; // discrepany between k-space resolution and image resolution - double yMaxFov = kyMax*upsampling; // actual FOV in y-direction (in x-direction xMax==FOV) + // k-space center at maximum echo + double fromMaxEcho = 0; + if ( (int)kyMax%2==0 ) + { + fromMaxEcho = -m_Parameters.m_SignalGen.m_tLine*(kyMax-1)/2 + dt*(kxMax-(int)kxMax%2)/2; + } + else + fromMaxEcho = -m_Parameters.m_SignalGen.m_tLine*(kyMax-1)/2 - dt*(kxMax-(int)kxMax%2)/2; - int xRingingOffset = xMax-kxMax; - int yRingingOffset = yMaxFov-kyMax; + double yMaxFov = yMax*m_Parameters.m_SignalGen.m_CroppingFactor; // actual FOV in y-direction (in x-direction FOV=xMax) - double noiseVar = m_Parameters.m_SignalGen.m_PartialFourier*m_Parameters.m_SignalGen.m_NoiseVariance/(yMaxFov*kxMax); // adjust noise variance since it is the intended variance in physical space and not in k-space + double noiseVar = m_Parameters.m_SignalGen.m_PartialFourier*m_Parameters.m_SignalGen.m_NoiseVariance/(kyMax*kxMax); // adjust noise variance since it is the intended variance in physical space and not in k-space while( !oit.IsAtEnd() ) { itk::Index< 2 > kIdx; kIdx[0] = oit.GetIndex()[0]; kIdx[1] = oit.GetIndex()[1]; // dephasing time double t= fromMaxEcho + ((double)kIdx[1]*kxMax+(double)kIdx[0])*dt; // readout time - double tall = 0; -// if (!m_Parameters.m_SignalGen.m_ReversePhase) - tall = ((double)kIdx[1]*kxMax+(double)kIdx[0])*dt; -// else -// tall = ((double)(kyMax-1-kIdx[1])*kxMax+(double)kIdx[0])*dt; + double tall = ((double)kIdx[1]*kxMax+(double)kIdx[0])*dt; // calculate eddy current decay factor double eddyDecay = 0; - if ( m_Parameters.m_SignalGen.m_EddyStrength>0) + if ( m_Parameters.m_Misc.m_CheckAddEddyCurrentsBox && m_Parameters.m_SignalGen.m_EddyStrength>0) eddyDecay = exp(-tall/m_Parameters.m_SignalGen.m_Tau ); // calcualte signal relaxation factors std::vector< double > relaxFactor; if ( m_Parameters.m_SignalGen.m_DoSimulateRelaxation) for (unsigned int i=0; ikyMax*m_Parameters.m_SignalGen.m_PartialFourier) pf = true; // reverse readout direction if (oit.GetIndex()[1]%2 == 1) kIdx[0] = kxMax-kIdx[0]-1; -// if (!pf) -// m_KSpaceImage->SetPixel(kIdx, t ); - - // rearrange slice - if( kIdx[0] < kxMax/2 ) - kIdx[0] = kIdx[0] + kxMax/2; + // shift k for DFT: (0 -- N) --> (-N/2 -- N/2) + double kx = kIdx[0]; + double ky = kIdx[1]; + if ((int)kxMax%2==1) + kx -= (kxMax-1)/2; else - kIdx[0] = kIdx[0] - kxMax/2; - - if( kIdx[1] < kyMax/2 ) - kIdx[1] = kIdx[1] + kyMax/2; + kx -= kxMax/2; + if ((int)kyMax%2==1) + ky -= (kyMax-1)/2; else - kIdx[1] = kIdx[1] - kyMax/2; + ky -= kyMax/2; // add ghosting - double kx = kIdx[0]; - double ky = kIdx[1]; if (oit.GetIndex()[1]%2 == 1) kx -= m_Parameters.m_SignalGen.m_KspaceLineOffset; // add gradient delay induced offset else kx += m_Parameters.m_SignalGen.m_KspaceLineOffset; // add gradient delay induced offset if (!pf) { - // add gibbs ringing offset (cropps k-space) - if (kx>=kxMax/2) - kx += xRingingOffset; - if (ky>=kyMax/2) - ky += yRingingOffset; - vcl_complex s(0,0); InputIteratorType it(m_CompartmentImages.at(0), m_CompartmentImages.at(0)->GetLargestPossibleRegion() ); while( !it.IsAtEnd() ) { - double x = it.GetIndex()[0]-xMax/2; - double y = it.GetIndex()[1]-yMax/2; + double x = it.GetIndex()[0]; + double y = it.GetIndex()[1]; + if ((int)xMax%2==1) + x -= (xMax-1)/2; + else + x -= xMax/2; + if ((int)yMax%2==1) + y -= (yMax-1)/2; + else + y -= yMax/2; DoubleVectorType pos; pos[0] = x; pos[1] = y; pos[2] = m_Z; pos = m_Transform*pos/1000; // vector from image center to current position (in meter) vcl_complex f(0, 0); // sum compartment signals and simulate relaxation for (unsigned int i=0; i( m_CompartmentImages.at(i)->GetPixel(it.GetIndex()) * relaxFactor.at(i) * m_Parameters.m_SignalGen.m_SignalScale, 0); else f += std::complex( m_CompartmentImages.at(i)->GetPixel(it.GetIndex()) * m_Parameters.m_SignalGen.m_SignalScale ); if (m_Parameters.m_SignalGen.m_CoilSensitivityProfile!=SignalGenerationParameters::COIL_CONSTANT) f *= CoilSensitivity(pos); // simulate eddy currents and other distortions double omega = 0; // frequency offset - if ( m_Parameters.m_SignalGen.m_EddyStrength>0 && !m_IsBaseline) + if ( m_Parameters.m_SignalGen.m_EddyStrength>0 && m_Parameters.m_Misc.m_CheckAddEddyCurrentsBox && !m_IsBaseline) { omega += (m_DiffusionGradientDirection[0]*pos[0]+m_DiffusionGradientDirection[1]*pos[1]+m_DiffusionGradientDirection[2]*pos[2]) * eddyDecay; } if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull()) // simulate distortions { itk::Point point3D; ItkDoubleImgType::IndexType index; index[0] = it.GetIndex()[0]; index[1] = it.GetIndex()[1]; index[2] = m_Zidx; if (m_Parameters.m_SignalGen.m_DoAddMotion) { m_Parameters.m_SignalGen.m_FrequencyMap->TransformIndexToPhysicalPoint(index, point3D); point3D = m_FiberBundle->TransformPoint(point3D.GetVnlVector(), -m_Rotation[0],-m_Rotation[1],-m_Rotation[2],-m_Translation[0],-m_Translation[1],-m_Translation[2]); m_Parameters.m_SignalGen.m_FrequencyMap->TransformPhysicalPointToIndex(point3D, index); if (m_Parameters.m_SignalGen.m_FrequencyMap->GetLargestPossibleRegion().IsInside(index)) omega += m_Parameters.m_SignalGen.m_FrequencyMap->GetPixel(index); } else { omega += m_Parameters.m_SignalGen.m_FrequencyMap->GetPixel(index); } } + // if signal comes from outside FOV, mirror it back (wrap-around artifact - aliasing) if (y<-yMaxFov/2) y += yMaxFov; else if (y>=yMaxFov/2) y -= yMaxFov; // actual DFT term s += f * exp( std::complex(0, 2 * M_PI * (kx*x/xMax + ky*y/yMaxFov + omega*t/1000 )) ); ++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_Parameters.m_SignalGen.m_NoiseVariance>0) s = vcl_complex(s.real()+randGen->GetNormalVariate(0,noiseVar), s.imag()+randGen->GetNormalVariate(0,noiseVar)); outputImage->SetPixel(kIdx, s); - m_KSpaceImage->SetPixel(oit.GetIndex(), sqrt(s.imag()*s.imag()+s.real()*s.real()) ); -// m_KSpaceImage->SetPixel(kIdx, sqrt(s.imag()*s.imag()+s.real()*s.real()) ); +// m_KSpaceImage->SetPixel(kIdx, t/dt ); + m_KSpaceImage->SetPixel(kIdx, sqrt(s.imag()*s.imag()+s.real()*s.real()) ); } ++oit; } } template< class TPixelType > void KspaceImageFilter< TPixelType > ::AfterThreadedGenerateData() { typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); double kxMax = outputImage->GetLargestPossibleRegion().GetSize(0); // k-space size in x-direction double kyMax = outputImage->GetLargestPossibleRegion().GetSize(1); // k-space size in y-direction ImageRegionIterator< OutputImageType > oit(outputImage, outputImage->GetLargestPossibleRegion()); while( !oit.IsAtEnd() ) // use hermitian k-space symmetry to fill empty k-space parts resulting from partial fourier acquisition { itk::Index< 2 > kIdx; kIdx[0] = oit.GetIndex()[0]; kIdx[1] = oit.GetIndex()[1]; // reverse phase if (!m_Parameters.m_SignalGen.m_ReversePhase) kIdx[1] = 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; - // flip k-space - if( kIdx[0] < kxMax/2 ) - kIdx[0] = kIdx[0] + kxMax/2; - else - kIdx[0] = kIdx[0] - kxMax/2; - - if( kIdx[1] < kyMax/2 ) - kIdx[1] = kIdx[1] + kyMax/2; - else - kIdx[1] = kIdx[1] - kyMax/2; - + // calculate symmetric index itk::Index< 2 > kIdx2; - kIdx2[0] = (int)(kxMax-kIdx[0])%(int)kxMax; - kIdx2[1] = (int)(kyMax-kIdx[1])%(int)kyMax; + kIdx2[0] = (int)(kxMax-kIdx[0]-(int)kxMax%2)%(int)kxMax; + kIdx2[1] = (int)(kyMax-kIdx[1]-(int)kyMax%2)%(int)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(oit.GetIndex(), sqrt(s.imag()*s.imag()+s.real()*s.real()) ); -// m_KSpaceImage->SetPixel(kIdx, sqrt(s.imag()*s.imag()+s.real()*s.real()) ); + 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()%(int)kxMax; spikeIdx[1] = randGen->GetIntegerVariate()%(int)kyMax; outputImage->SetPixel(spikeIdx, m_Spike); } } } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp index 83fabc3d58..c0532267c6 100755 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp @@ -1,1568 +1,1569 @@ /*=================================================================== 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 namespace itk { template< class PixelType > TractsToDWIImageFilter< PixelType >::TractsToDWIImageFilter() : m_FiberBundle(NULL) , m_StatusText("") , m_UseConstantRandSeed(false) , m_RandGen(itk::Statistics::MersenneTwisterRandomVariateGenerator::New()) { m_RandGen->SetSeed(); } template< class PixelType > TractsToDWIImageFilter< PixelType >::~TractsToDWIImageFilter() { } template< class PixelType > TractsToDWIImageFilter< PixelType >::DoubleDwiType::Pointer TractsToDWIImageFilter< PixelType >::SimulateKspaceAcquisition( std::vector< DoubleDwiType::Pointer >& images ) { int numFiberCompartments = m_Parameters.m_FiberModelList.size(); // create slice object ImageRegion<2> sliceRegion; sliceRegion.SetSize(0, m_WorkingImageRegion.GetSize()[0]); sliceRegion.SetSize(1, m_WorkingImageRegion.GetSize()[1]); Vector< double, 2 > sliceSpacing; sliceSpacing[0] = m_WorkingSpacing[0]; sliceSpacing[1] = m_WorkingSpacing[1]; DoubleDwiType::PixelType nullPix; nullPix.SetSize(images.at(0)->GetVectorLength()); nullPix.Fill(0.0); DoubleDwiType::Pointer 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_CroppedRegion ); magnitudeDwiImage->SetBufferedRegion( m_CroppedRegion ); magnitudeDwiImage->SetRequestedRegion( m_CroppedRegion ); magnitudeDwiImage->SetVectorLength( 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_CroppedRegion ); m_PhaseImage->SetBufferedRegion( m_CroppedRegion ); m_PhaseImage->SetRequestedRegion( m_CroppedRegion ); m_PhaseImage->SetVectorLength( images.at(0)->GetVectorLength() ); m_PhaseImage->Allocate(); m_PhaseImage->FillBuffer(nullPix); // itk::ImageRegion<4> imageRegion4D; // itk::Vector imageSpacing4D; imageSpacing4D.Fill(1); // itk::Point imageOrigin4D; imageOrigin4D.Fill(0); // itk::Matrix imageDirection4D; imageDirection4D.SetIdentity(); // imageRegion4D.SetSize(0, m_CroppedRegion.GetSize(0)); // imageRegion4D.SetSize(1, m_CroppedRegion.GetSize(1)); // imageRegion4D.SetSize(2, m_CroppedRegion.GetSize(2)); // imageRegion4D.SetSize(3, m_Parameters.m_SignalGen.m_NumberOfCoils); // for (int i=0; i<3; i++) // { // imageSpacing4D[i] = m_Parameters.m_SignalGen.m_ImageSpacing[i]; // imageOrigin4D[i] = m_Parameters.m_SignalGen.m_ImageOrigin[i]; // for (int j=0; j<3; j++) // imageDirection4D[i][j]=m_Parameters.m_SignalGen.m_ImageDirection[i][j]; // } //// ItkDoubleImgType4D::PixelType nullPix4D; //// nullPix4D.SetSize(images.at(0)->GetVectorLength()); //// nullPix4D.Fill(0.0); // m_KspaceImage = ItkDoubleImgType4D::New(); // m_KspaceImage->SetSpacing( imageSpacing4D ); // m_KspaceImage->SetOrigin( imageOrigin4D ); // m_KspaceImage->SetDirection( imageDirection4D ); // m_KspaceImage->SetLargestPossibleRegion( imageRegion4D ); // m_KspaceImage->SetBufferedRegion( imageRegion4D ); // m_KspaceImage->SetRequestedRegion( imageRegion4D ); //// m_KspaceImage->SetVectorLength( images.at(0)->GetVectorLength() ); // m_KspaceImage->Allocate(); // m_KspaceImage->FillBuffer(0); m_KspaceImage = DoubleDwiType::New(); m_KspaceImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); m_KspaceImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin ); m_KspaceImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); m_KspaceImage->SetLargestPossibleRegion( m_CroppedRegion ); m_KspaceImage->SetBufferedRegion( m_CroppedRegion ); m_KspaceImage->SetRequestedRegion( m_CroppedRegion ); m_KspaceImage->SetVectorLength( m_Parameters.m_SignalGen.m_NumberOfCoils ); m_KspaceImage->Allocate(); m_KspaceImage->FillBuffer(nullPix); std::vector< unsigned int > spikeVolume; for (unsigned int i=0; iGetIntegerVariate()%(images.at(0)->GetVectorLength())); std::sort (spikeVolume.begin(), spikeVolume.end()); std::reverse (spikeVolume.begin(), spikeVolume.end()); // 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 * M_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()); } m_StatusText += "0% 10 20 30 40 50 60 70 80 90 100%\n"; m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*"; unsigned long lastTick = 0; boost::progress_display disp(images.at(0)->GetVectorLength()*images.at(0)->GetLargestPossibleRegion().GetSize(2)); for (unsigned int g=0; gGetVectorLength(); g++) { std::vector< unsigned int > spikeSlice; while (!spikeVolume.empty() && spikeVolume.back()==g) { spikeSlice.push_back(m_RandGen->GetIntegerVariate()%images.at(0)->GetLargestPossibleRegion().GetSize(2)); spikeVolume.pop_back(); } std::sort (spikeSlice.begin(), spikeSlice.end()); std::reverse (spikeSlice.begin(), spikeSlice.end()); for (unsigned int z=0; zGetLargestPossibleRegion().GetSize(2); z++) { std::vector< SliceType::Pointer > compartmentSlices; std::vector< double > t2Vector; std::vector< double > t1Vector; for (unsigned int i=0; i* signalModel; if (iSetLargestPossibleRegion( sliceRegion ); slice->SetBufferedRegion( sliceRegion ); slice->SetRequestedRegion( sliceRegion ); slice->SetSpacing(sliceSpacing); slice->Allocate(); slice->FillBuffer(0.0); // extract slice from channel g for (unsigned int y=0; yGetLargestPossibleRegion().GetSize(1); y++) for (unsigned int x=0; xGetLargestPossibleRegion().GetSize(0); x++) { SliceType::IndexType index2D; index2D[0]=x; index2D[1]=y; DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z; slice->SetPixel(index2D, images.at(i)->GetPixel(index3D)[g]); } compartmentSlices.push_back(slice); t2Vector.push_back(signalModel->GetT2()); t1Vector.push_back(signalModel->GetT1()); } int numSpikes = 0; while (!spikeSlice.empty() && spikeSlice.back()==z) { numSpikes++; spikeSlice.pop_back(); } int spikeCoil = m_RandGen->GetIntegerVariate()%m_Parameters.m_SignalGen.m_NumberOfCoils; if (this->GetAbortGenerateData()) return NULL; #pragma omp parallel for for (int c=0; c outSize; outSize.SetElement(0, m_CroppedRegion.GetSize(0)); outSize.SetElement(1, m_CroppedRegion.GetSize(1)); itk::KspaceImageFilter< SliceType::PixelType >::Pointer idft = itk::KspaceImageFilter< SliceType::PixelType >::New(); idft->SetCompartmentImages(compartmentSlices); idft->SetT2(t2Vector); idft->SetT1(t1Vector); idft->SetUseConstantRandSeed(m_UseConstantRandSeed); idft->SetParameters(m_Parameters); idft->SetZ((double)z-(double)(images.at(0)->GetLargestPossibleRegion().GetSize(2)-images.at(0)->GetLargestPossibleRegion().GetSize(2)%2)/2.0); idft->SetZidx(z); idft->SetCoilPosition(coilPositions.at(c)); idft->SetFiberBundle(m_FiberBundleWorkingCopy); if (m_Parameters.m_SignalGen.m_DoAddMotion) { idft->SetTranslation(m_Translations.at(g)); idft->SetRotation(m_Rotations.at(g)); } idft->SetDiffusionGradientDirection(m_Parameters.m_SignalGen.GetGradientDirection(g)); idft->SetOutSize(outSize); if (c==spikeCoil) idft->SetSpikesPerSlice(numSpikes); // idft->SetNumberOfThreads(1); idft->Update(); ComplexSliceType::Pointer fSlice; fSlice = idft->GetOutput(); // fourier transform slice ComplexSliceType::Pointer newSlice; itk::DftImageFilter< SliceType::PixelType >::Pointer dft = itk::DftImageFilter< SliceType::PixelType >::New(); dft->SetInput(fSlice); dft->SetParameters(m_Parameters); dft->Update(); newSlice = dft->GetOutput(); // put slice back into channel g for (unsigned int y=0; 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; ComplexSliceType::IndexType index2D; index2D[0]=x; index2D[1]=y; ComplexSliceType::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 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 // { // DoubleDwiType4D::IndexType idx4d; // idx4d[0]=index3D[0]; // idx4d[1]=index3D[1]; // idx4d[2]=index3D[2]; // idx4d[3]=c; // ItkDoubleImgType4D::PixelType pix4D = m_KspaceImage->GetPixel(idx4d); // pix4D = idft->GetKSpaceImage()->GetPixel(index2D); // m_KspaceImage->SetPixel(idx4d, pix4D); // } 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) { #pragma omp parallel for collapse(2) 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; 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++) m_StatusText += "*"; lastTick = newTick; } } m_StatusText += "\n\n"; 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()>900) // it.Set(900); if (it.Get()>max) max = it.Get(); if (it.Get()::Pointer scaler = itk::ShiftScaleImageFilter< ItkDoubleImgType, ItkDoubleImgType >::New(); scaler->SetInput(image); scaler->SetShift(-min); scaler->SetScale(1.0/(max-min)); scaler->Update(); return scaler->GetOutput(); } template< class PixelType > void TractsToDWIImageFilter< PixelType >::CheckVolumeFractionImages() { m_UseRelativeNonFiberVolumeFractions = false; // check for fiber volume fraction maps int fibVolImages = 0; for (int i=0; iGetVolumeFractionImage().IsNotNull()) { m_StatusText += "Using volume fraction map for fiber compartment " + boost::lexical_cast(i+1) + "\n"; MITK_INFO << "Using volume fraction map for fiber compartment " + boost::lexical_cast(i+1); fibVolImages++; } // check for non-fiber volume fraction maps int nonfibVolImages = 0; for (int i=0; iGetVolumeFractionImage().IsNotNull()) { m_StatusText += "Using volume fraction map for non-fiber compartment " + boost::lexical_cast(i+1) + "\n"; MITK_INFO << "Using volume fraction map for non-fiber compartment " + boost::lexical_cast(i+1); nonfibVolImages++; } // not all fiber compartments are using volume fraction maps --> non-fiber volume fractions are assumed to be relative to the non-fiber volume and not absolute voxel-volume fractions. // this means if two non-fiber compartments are used but only one of them has an associated volume fraction map, the repesctive other volume fraction map can be determined as inverse (1-val) of the present volume fraction map- if ( fibVolImages::Pointer inverter = itk::InvertIntensityImageFilter< ItkDoubleImgType, ItkDoubleImgType >::New(); inverter->SetMaximum(1.0); if ( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage().IsNull() && m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage().IsNotNull() ) { // m_Parameters.m_NonFiberModelList[1]->SetVolumeFractionImage( NormalizeInsideMask( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage() ) ); inverter->SetInput( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage() ); inverter->Update(); m_Parameters.m_NonFiberModelList[0]->SetVolumeFractionImage(inverter->GetOutput()); } else if ( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage().IsNull() && m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage().IsNotNull() ) { // m_Parameters.m_NonFiberModelList[0]->SetVolumeFractionImage( NormalizeInsideMask( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage() ) ); inverter->SetInput( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage() ); inverter->Update(); m_Parameters.m_NonFiberModelList[1]->SetVolumeFractionImage(inverter->GetOutput()); } else { itkExceptionMacro("Something went wrong in automatically calculating the missing non-fiber volume fraction image! Did you use two non fiber compartments but only one volume fraction image? Then it should work and this error is really strange."); } 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) { m_StatusText += "Not all fiber compartments are using an associated volume fraction image.\nAssuming non-fiber volume fraction images to contain values relative to the remaining non-fiber volume, not absolute values.\n"; MITK_INFO << "Not all fiber compartments are using an associated volume fraction image.\nAssuming non-fiber volume fraction images to contain values relative to the remaining non-fiber volume, not absolute values."; m_UseRelativeNonFiberVolumeFractions = true; // 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 (int i=0; iSetSpacing( m_WorkingSpacing ); doubleImg->SetOrigin( m_WorkingOrigin ); doubleImg->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); doubleImg->SetLargestPossibleRegion( m_WorkingImageRegion ); doubleImg->SetBufferedRegion( m_WorkingImageRegion ); doubleImg->SetRequestedRegion( m_WorkingImageRegion ); doubleImg->Allocate(); doubleImg->FillBuffer(0); m_VolumeFractions.push_back(doubleImg); } } template< class PixelType > void TractsToDWIImageFilter< PixelType >::InitializeData() { // initialize output dwi image m_CroppedRegion = m_Parameters.m_SignalGen.m_ImageRegion; m_CroppedRegion.SetSize(1, m_CroppedRegion.GetSize(1)*m_Parameters.m_SignalGen.m_CroppingFactor); itk::Point shiftedOrigin = m_Parameters.m_SignalGen.m_ImageOrigin; shiftedOrigin[1] += (m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1)-m_CroppedRegion.GetSize(1))*m_Parameters.m_SignalGen.m_ImageSpacing[1]/2; m_OutputImage = OutputImageType::New(); m_OutputImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); m_OutputImage->SetOrigin( shiftedOrigin ); m_OutputImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); m_OutputImage->SetLargestPossibleRegion( m_CroppedRegion ); m_OutputImage->SetBufferedRegion( m_CroppedRegion ); m_OutputImage->SetRequestedRegion( m_CroppedRegion ); m_OutputImage->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() ); m_OutputImage->Allocate(); typename OutputImageType::PixelType temp; temp.SetSize(m_Parameters.m_SignalGen.GetNumVolumes()); temp.Fill(0.0); m_OutputImage->FillBuffer(temp); - if ( m_CroppedRegion.GetSize(0)%2 == 1 ) - m_CroppedRegion.SetSize(0, m_CroppedRegion.GetSize(0)+1); - if ( m_CroppedRegion.GetSize(1)%2 == 1 ) - m_CroppedRegion.SetSize(1, m_CroppedRegion.GetSize(1)+1); - - // ADJUST GEOMETRY FOR FURTHER PROCESSING - // is input slize size a power of two? - unsigned int x=m_Parameters.m_SignalGen.m_ImageRegion.GetSize(0); unsigned int y=m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1); - ItkDoubleImgType::SizeType pad; pad[0]=x%2; pad[1]=y%2; pad[2]=0; - m_Parameters.m_SignalGen.m_ImageRegion.SetSize(0, x+pad[0]); - m_Parameters.m_SignalGen.m_ImageRegion.SetSize(1, y+pad[1]); - if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull() && (pad[0]>0 || pad[1]>0)) - { - itk::ConstantPadImageFilter::Pointer zeroPadder = itk::ConstantPadImageFilter::New(); - zeroPadder->SetInput(m_Parameters.m_SignalGen.m_FrequencyMap); - zeroPadder->SetConstant(0); - zeroPadder->SetPadUpperBound(pad); - zeroPadder->Update(); - m_Parameters.m_SignalGen.m_FrequencyMap = zeroPadder->GetOutput(); - } - if (m_Parameters.m_SignalGen.m_MaskImage.IsNotNull() && (pad[0]>0 || pad[1]>0)) - { - itk::ConstantPadImageFilter::Pointer zeroPadder = itk::ConstantPadImageFilter::New(); - zeroPadder->SetInput(m_Parameters.m_SignalGen.m_MaskImage); - zeroPadder->SetConstant(0); - zeroPadder->SetPadUpperBound(pad); - zeroPadder->Update(); - m_Parameters.m_SignalGen.m_MaskImage = zeroPadder->GetOutput(); - } + //// NOT NECESSARY ANYMORE +// if ( m_CroppedRegion.GetSize(0)%2 == 1 ) +// m_CroppedRegion.SetSize(0, m_CroppedRegion.GetSize(0)+1); +// if ( m_CroppedRegion.GetSize(1)%2 == 1 ) +// m_CroppedRegion.SetSize(1, m_CroppedRegion.GetSize(1)+1); + +// // ADJUST GEOMETRY FOR FURTHER PROCESSING +// // is input slize size a power of two? +// unsigned int x=m_Parameters.m_SignalGen.m_ImageRegion.GetSize(0); unsigned int y=m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1); +// ItkDoubleImgType::SizeType pad; pad[0]=x%2; pad[1]=y%2; pad[2]=0; +// m_Parameters.m_SignalGen.m_ImageRegion.SetSize(0, x+pad[0]); +// m_Parameters.m_SignalGen.m_ImageRegion.SetSize(1, y+pad[1]); +// if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull() && (pad[0]>0 || pad[1]>0)) +// { +// itk::ConstantPadImageFilter::Pointer zeroPadder = itk::ConstantPadImageFilter::New(); +// zeroPadder->SetInput(m_Parameters.m_SignalGen.m_FrequencyMap); +// zeroPadder->SetConstant(0); +// zeroPadder->SetPadUpperBound(pad); +// zeroPadder->Update(); +// m_Parameters.m_SignalGen.m_FrequencyMap = zeroPadder->GetOutput(); +// } +// if (m_Parameters.m_SignalGen.m_MaskImage.IsNotNull() && (pad[0]>0 || pad[1]>0)) +// { +// itk::ConstantPadImageFilter::Pointer zeroPadder = itk::ConstantPadImageFilter::New(); +// zeroPadder->SetInput(m_Parameters.m_SignalGen.m_MaskImage); +// zeroPadder->SetConstant(0); +// zeroPadder->SetPadUpperBound(pad); +// zeroPadder->Update(); +// m_Parameters.m_SignalGen.m_MaskImage = zeroPadder->GetOutput(); +// } // Apply in-plane upsampling for Gibbs ringing artifact double upsampling = 1; if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging) upsampling = 2; m_WorkingSpacing = m_Parameters.m_SignalGen.m_ImageSpacing; m_WorkingSpacing[0] /= upsampling; m_WorkingSpacing[1] /= upsampling; m_WorkingImageRegion = m_Parameters.m_SignalGen.m_ImageRegion; m_WorkingImageRegion.SetSize(0, m_Parameters.m_SignalGen.m_ImageRegion.GetSize()[0]*upsampling); m_WorkingImageRegion.SetSize(1, m_Parameters.m_SignalGen.m_ImageRegion.GetSize()[1]*upsampling); m_WorkingOrigin = m_Parameters.m_SignalGen.m_ImageOrigin; m_WorkingOrigin[0] -= m_Parameters.m_SignalGen.m_ImageSpacing[0]/2; m_WorkingOrigin[0] += m_WorkingSpacing[0]/2; m_WorkingOrigin[1] -= m_Parameters.m_SignalGen.m_ImageSpacing[1]/2; m_WorkingOrigin[1] += m_WorkingSpacing[1]/2; m_WorkingOrigin[2] -= m_Parameters.m_SignalGen.m_ImageSpacing[2]/2; m_WorkingOrigin[2] += m_WorkingSpacing[2]/2; m_VoxelVolume = m_WorkingSpacing[0]*m_WorkingSpacing[1]*m_WorkingSpacing[2]; // generate double images to store the individual compartment signals m_CompartmentImages.clear(); int numFiberCompartments = m_Parameters.m_FiberModelList.size(); int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size(); for (int i=0; iSetSpacing( m_WorkingSpacing ); doubleDwi->SetOrigin( m_WorkingOrigin ); doubleDwi->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); doubleDwi->SetLargestPossibleRegion( m_WorkingImageRegion ); doubleDwi->SetBufferedRegion( m_WorkingImageRegion ); doubleDwi->SetRequestedRegion( m_WorkingImageRegion ); doubleDwi->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() ); doubleDwi->Allocate(); DoubleDwiType::PixelType pix; pix.SetSize(m_Parameters.m_SignalGen.GetNumVolumes()); pix.Fill(0.0); doubleDwi->FillBuffer(pix); m_CompartmentImages.push_back(doubleDwi); } // resample mask image and frequency map to fit upsampled geometry if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging) { if (m_Parameters.m_SignalGen.m_MaskImage.IsNotNull()) { // rescale mask image (otherwise there are problems with the resampling) itk::RescaleIntensityImageFilter::Pointer rescaler = itk::RescaleIntensityImageFilter::New(); rescaler->SetInput(0,m_Parameters.m_SignalGen.m_MaskImage); rescaler->SetOutputMaximum(100); rescaler->SetOutputMinimum(0); rescaler->Update(); // resample mask image itk::ResampleImageFilter::Pointer resampler = itk::ResampleImageFilter::New(); resampler->SetInput(rescaler->GetOutput()); resampler->SetOutputParametersFromImage(m_Parameters.m_SignalGen.m_MaskImage); resampler->SetSize(m_WorkingImageRegion.GetSize()); resampler->SetOutputSpacing(m_WorkingSpacing); resampler->SetOutputOrigin(m_WorkingOrigin); itk::NearestNeighborInterpolateImageFunction::Pointer nn_interpolator = itk::NearestNeighborInterpolateImageFunction::New(); resampler->SetInterpolator(nn_interpolator); resampler->Update(); m_Parameters.m_SignalGen.m_MaskImage = resampler->GetOutput(); } // resample frequency map if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull()) { itk::ResampleImageFilter::Pointer resampler = itk::ResampleImageFilter::New(); resampler->SetInput(m_Parameters.m_SignalGen.m_FrequencyMap); resampler->SetOutputParametersFromImage(m_Parameters.m_SignalGen.m_FrequencyMap); resampler->SetSize(m_WorkingImageRegion.GetSize()); resampler->SetOutputSpacing(m_WorkingSpacing); resampler->SetOutputOrigin(m_WorkingOrigin); itk::NearestNeighborInterpolateImageFunction::Pointer nn_interpolator = itk::NearestNeighborInterpolateImageFunction::New(); resampler->SetInterpolator(nn_interpolator); resampler->Update(); m_Parameters.m_SignalGen.m_FrequencyMap = resampler->GetOutput(); } } m_MaskImageSet = true; if (m_Parameters.m_SignalGen.m_MaskImage.IsNull()) { // no input tissue mask is set -> create default m_StatusText += "No tissue mask set\n"; MITK_INFO << "No tissue mask set"; m_Parameters.m_SignalGen.m_MaskImage = ItkUcharImgType::New(); m_Parameters.m_SignalGen.m_MaskImage->SetSpacing( m_WorkingSpacing ); m_Parameters.m_SignalGen.m_MaskImage->SetOrigin( m_WorkingOrigin ); m_Parameters.m_SignalGen.m_MaskImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); m_Parameters.m_SignalGen.m_MaskImage->SetLargestPossibleRegion( m_WorkingImageRegion ); m_Parameters.m_SignalGen.m_MaskImage->SetBufferedRegion( m_WorkingImageRegion ); m_Parameters.m_SignalGen.m_MaskImage->SetRequestedRegion( m_WorkingImageRegion ); m_Parameters.m_SignalGen.m_MaskImage->Allocate(); m_Parameters.m_SignalGen.m_MaskImage->FillBuffer(100); m_MaskImageSet = false; } else { if (m_Parameters.m_SignalGen.m_MaskImage->GetLargestPossibleRegion()!=m_WorkingImageRegion) itkExceptionMacro("Mask image and specified DWI geometry are not matching!"); m_StatusText += "Using tissue mask\n"; MITK_INFO << "Using tissue mask"; } if (m_Parameters.m_SignalGen.m_DoAddMotion) { std::string fileName = "fiberfox_motion_0.log"; std::string filePath = mitk::IOUtil::GetTempPath(); if (m_Parameters.m_Misc.m_OutputPath.size()>0) filePath = m_Parameters.m_Misc.m_OutputPath; int c = 1; while (itksys::SystemTools::FileExists((filePath+fileName).c_str())) { fileName = "fiberfox_motion_"; fileName += boost::lexical_cast(c); fileName += ".log"; c++; } m_Logfile.open((filePath+fileName).c_str()); m_Logfile << "0 rotation: 0,0,0; translation: 0,0,0\n"; if (m_Parameters.m_SignalGen.m_DoRandomizeMotion) { m_StatusText += "Adding random motion artifacts:\n"; m_StatusText += "Maximum rotation: +/-" + boost::lexical_cast(m_Parameters.m_SignalGen.m_Rotation) + "°\n"; m_StatusText += "Maximum translation: +/-" + boost::lexical_cast(m_Parameters.m_SignalGen.m_Translation) + "mm\n"; } else { m_StatusText += "Adding linear motion artifacts:\n"; m_StatusText += "Maximum rotation: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_Rotation) + "°\n"; m_StatusText += "Maximum translation: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_Translation) + "mm\n"; } m_StatusText += "Motion logfile: " + (filePath+fileName) + "\n"; MITK_INFO << "Adding motion artifacts"; MITK_INFO << "Maximum rotation: " << m_Parameters.m_SignalGen.m_Rotation; MITK_INFO << "Maxmimum translation: " << m_Parameters.m_SignalGen.m_Translation; MITK_INFO << "Motion logfile: " << filePath << fileName; } // get transform for motion artifacts m_Rotation.Fill(0.0); m_Translation.Fill(0.0); m_Rotations.push_back(m_Rotation); m_Translations.push_back(m_Translation); m_Rotation = m_Parameters.m_SignalGen.m_Rotation/m_Parameters.m_SignalGen.GetNumVolumes(); m_Translation = m_Parameters.m_SignalGen.m_Translation/m_Parameters.m_SignalGen.GetNumVolumes(); // creat image to hold transformed mask (motion artifact) m_TransformedMaskImage = ItkUcharImgType::New(); itk::ImageDuplicator::Pointer duplicator = itk::ImageDuplicator::New(); duplicator->SetInputImage(m_Parameters.m_SignalGen.m_MaskImage); duplicator->Update(); m_TransformedMaskImage = duplicator->GetOutput(); // second upsampling needed for motion artifacts ImageRegion<3> upsampledImageRegion = m_WorkingImageRegion; DoubleVectorType upsampledSpacing = m_WorkingSpacing; upsampledSpacing[0] /= 4; upsampledSpacing[1] /= 4; upsampledSpacing[2] /= 4; upsampledImageRegion.SetSize(0, m_WorkingImageRegion.GetSize()[0]*4); upsampledImageRegion.SetSize(1, m_WorkingImageRegion.GetSize()[1]*4); upsampledImageRegion.SetSize(2, m_WorkingImageRegion.GetSize()[2]*4); itk::Point upsampledOrigin = m_WorkingOrigin; upsampledOrigin[0] -= m_WorkingSpacing[0]/2; upsampledOrigin[0] += upsampledSpacing[0]/2; upsampledOrigin[1] -= m_WorkingSpacing[1]/2; upsampledOrigin[1] += upsampledSpacing[1]/2; upsampledOrigin[2] -= m_WorkingSpacing[2]/2; upsampledOrigin[2] += upsampledSpacing[2]/2; m_UpsampledMaskImage = ItkUcharImgType::New(); itk::ResampleImageFilter::Pointer upsampler = itk::ResampleImageFilter::New(); upsampler->SetInput(m_Parameters.m_SignalGen.m_MaskImage); upsampler->SetOutputParametersFromImage(m_Parameters.m_SignalGen.m_MaskImage); upsampler->SetSize(upsampledImageRegion.GetSize()); upsampler->SetOutputSpacing(upsampledSpacing); upsampler->SetOutputOrigin(upsampledOrigin); itk::NearestNeighborInterpolateImageFunction::Pointer nn_interpolator = itk::NearestNeighborInterpolateImageFunction::New(); upsampler->SetInterpolator(nn_interpolator); upsampler->Update(); m_UpsampledMaskImage = upsampler->GetOutput(); } template< class PixelType > void TractsToDWIImageFilter< PixelType >::InitializeFiberData() { // resample fiber bundle for sufficient voxel coverage m_StatusText += "\n"+this->GetTime()+" > Resampling fibers ...\n"; m_SegmentVolume = 0.0001; float minSpacing = 1; if(m_WorkingSpacing[0]GetDeepCopy(); // working copy is needed because we need to resample the fibers but do not want to change the original bundle double volumeAccuracy = 10; m_FiberBundleWorkingCopy->ResampleSpline(minSpacing/volumeAccuracy); m_mmRadius = m_Parameters.m_SignalGen.m_AxonRadius/1000; if (m_mmRadius>0) m_SegmentVolume = M_PI*m_mmRadius*m_mmRadius*minSpacing/volumeAccuracy; m_FiberBundleTransformed = m_FiberBundleWorkingCopy; // a secon fiber bundle is needed to store the transformed version of the m_FiberBundleWorkingCopy } template< class PixelType > void TractsToDWIImageFilter< PixelType >::GenerateData() { m_TimeProbe.Start(); m_StatusText = "Starting simulation\n"; // check input data if (m_FiberBundle.IsNull()) itkExceptionMacro("Input fiber bundle is NULL!"); if (m_Parameters.m_FiberModelList.empty()) itkExceptionMacro("No diffusion model for fiber compartments defined! At least one fiber compartment is necessary."); if (m_Parameters.m_NonFiberModelList.empty()) itkExceptionMacro("No diffusion model for non-fiber compartments defined! At least one non-fiber compartment is necessary."); int baselineIndex = m_Parameters.m_SignalGen.GetFirstBaselineIndex(); if (baselineIndex<0) itkExceptionMacro("No baseline index found!"); if (!m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition) // No upsampling of input image needed if no k-space simulation is performed m_Parameters.m_SignalGen.m_DoAddGibbsRinging = false; if (m_UseConstantRandSeed) // always generate the same random numbers? m_RandGen->SetSeed(0); else m_RandGen->SetSeed(); InitializeData(); CheckVolumeFractionImages(); InitializeFiberData(); int numFiberCompartments = m_Parameters.m_FiberModelList.size(); int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size(); double maxVolume = 0; unsigned long lastTick = 0; int signalModelSeed = m_RandGen->GetIntegerVariate(); m_StatusText += "\n"+this->GetTime()+" > Generating " + boost::lexical_cast(numFiberCompartments+numNonFiberCompartments) + "-compartment diffusion-weighted signal.\n"; MITK_INFO << "Generating " << numFiberCompartments+numNonFiberCompartments << "-compartment diffusion-weighted signal."; int numFibers = m_FiberBundleWorkingCopy->GetNumFibers(); boost::progress_display disp(numFibers*m_Parameters.m_SignalGen.GetNumVolumes()); switch (m_Parameters.m_SignalGen.m_DiffusionDirectionMode) { case(SignalGenerationParameters::FIBER_TANGENT_DIRECTIONS): // use fiber tangent directions to determine diffusion direction { m_StatusText += "0% 10 20 30 40 50 60 70 80 90 100%\n"; m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*"; for (unsigned int g=0; gSetSeed(signalModelSeed); for (int i=0; iSetSeed(signalModelSeed); // storing voxel-wise intra-axonal volume in mm³ ItkDoubleImgType::Pointer intraAxonalVolumeImage = ItkDoubleImgType::New(); intraAxonalVolumeImage->SetSpacing( m_WorkingSpacing ); intraAxonalVolumeImage->SetOrigin( m_WorkingOrigin ); intraAxonalVolumeImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); intraAxonalVolumeImage->SetLargestPossibleRegion( m_WorkingImageRegion ); intraAxonalVolumeImage->SetBufferedRegion( m_WorkingImageRegion ); intraAxonalVolumeImage->SetRequestedRegion( m_WorkingImageRegion ); intraAxonalVolumeImage->Allocate(); intraAxonalVolumeImage->FillBuffer(0); maxVolume = 0; vtkPolyData* fiberPolyData = m_FiberBundleTransformed->GetFiberPolyData(); // generate fiber signal (if there are any fiber models present) if (!m_Parameters.m_FiberModelList.empty()) for( int i=0; iGetFiberWeight(i); vtkCell* cell = fiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (numPoints<2) continue; for( int j=0; jGetAbortGenerateData()) { m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n"; return; } double* temp = points->GetPoint(j); itk::Point vertex = GetItkPoint(temp); itk::Vector v = GetItkVector(temp); itk::Vector dir(3); if (jGetPoint(j+1))-v; else dir = v-GetItkVector(points->GetPoint(j-1)); if (dir.GetSquaredNorm()<0.0001 || dir[0]!=dir[0] || dir[1]!=dir[1] || dir[2]!=dir[2]) continue; itk::Index<3> idx; itk::ContinuousIndex contIndex; m_TransformedMaskImage->TransformPhysicalPointToIndex(vertex, idx); m_TransformedMaskImage->TransformPhysicalPointToContinuousIndex(vertex, contIndex); if (!m_TransformedMaskImage->GetLargestPossibleRegion().IsInside(idx) || m_TransformedMaskImage->GetPixel(idx)<=0) continue; // generate signal for each fiber compartment for (int k=0; kSetFiberDirection(dir); DoubleDwiType::PixelType pix = m_CompartmentImages.at(k)->GetPixel(idx); pix[g] += fiberWeight*m_SegmentVolume*m_Parameters.m_FiberModelList[k]->SimulateMeasurement(g); m_CompartmentImages.at(k)->SetPixel(idx, pix); } // update fiber volume image double vol = intraAxonalVolumeImage->GetPixel(idx) + m_SegmentVolume*fiberWeight; intraAxonalVolumeImage->SetPixel(idx, vol); if (vol>maxVolume) // we assume that the first volume is always unweighted! maxVolume = vol; } // progress report ++disp; unsigned long newTick = 50*disp.count()/disp.expected_count(); for (unsigned int tick = 0; tick<(newTick-lastTick); tick++) m_StatusText += "*"; lastTick = newTick; } // generate non-fiber signal ImageRegionIterator it3(m_TransformedMaskImage, m_TransformedMaskImage->GetLargestPossibleRegion()); double fact = 1; // density correction factor in mm³ if (m_Parameters.m_SignalGen.m_AxonRadius<0.0001 || maxVolume>m_VoxelVolume) // the fullest voxel is always completely full fact = m_VoxelVolume/maxVolume; while(!it3.IsAtEnd()) { if (it3.Get()>0) { DoubleDwiType::IndexType index = it3.GetIndex(); itk::Point point; m_TransformedMaskImage->TransformIndexToPhysicalPoint(index, point); if (m_Parameters.m_SignalGen.m_DoAddMotion) { if (m_Parameters.m_SignalGen.m_DoRandomizeMotion && g>0) point = m_FiberBundleWorkingCopy->TransformPoint(point.GetVnlVector(), -m_Rotation[0],-m_Rotation[1],-m_Rotation[2],-m_Translation[0],-m_Translation[1],-m_Translation[2]); else if (g>=0) point = m_FiberBundleWorkingCopy->TransformPoint(point.GetVnlVector(), -m_Rotation[0]*g,-m_Rotation[1]*g,-m_Rotation[2]*g,-m_Translation[0]*g,-m_Translation[1]*g,-m_Translation[2]*g); } double iAxVolume = intraAxonalVolumeImage->GetPixel(index); // if volume fraction image is set use it, otherwise use scaling factor to obtain one full fiber voxel double fact2 = fact; if (m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()!=nullptr && iAxVolume>0.0001) { DoubleDwiType::IndexType newIndex; m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()->TransformPhysicalPointToIndex(point, newIndex); if (m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()->GetLargestPossibleRegion().IsInside(newIndex)) fact2 = m_VoxelVolume*m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()->GetPixel(newIndex)/iAxVolume; } // adjust intra-axonal image value for (int i=0; iGetPixel(index); pix[g] *= fact2; m_CompartmentImages.at(i)->SetPixel(index, pix); } // simulate other compartments SimulateExtraAxonalSignal(index, iAxVolume*fact2, g); } ++it3; } // move fibers SimulateMotion(g); } break; } case (SignalGenerationParameters::MAIN_FIBER_DIRECTIONS): // use main fiber directions to determine voxel-wise diffusion directions { typedef itk::Image< itk::Vector< float, 3>, 3 > ItkDirectionImage3DType; typedef itk::VectorContainer< unsigned int, ItkDirectionImage3DType::Pointer > ItkDirectionImageContainerType; // calculate main fiber directions itk::TractsToVectorImageFilter::Pointer fOdfFilter = itk::TractsToVectorImageFilter::New(); fOdfFilter->SetFiberBundle(m_FiberBundleTransformed); fOdfFilter->SetMaskImage(m_TransformedMaskImage); fOdfFilter->SetAngularThreshold(cos(m_Parameters.m_SignalGen.m_FiberSeparationThreshold*M_PI/180.0)); fOdfFilter->SetNormalizeVectors(false); fOdfFilter->SetUseWorkingCopy(true); fOdfFilter->SetSizeThreshold(0); fOdfFilter->SetMaxNumDirections(3); fOdfFilter->Update(); ItkDirectionImageContainerType::Pointer directionImageContainer = fOdfFilter->GetDirectionImageContainer(); // allocate image storing intra-axonal volume fraction information ItkDoubleImgType::Pointer intraAxonalVolumeImage = ItkDoubleImgType::New(); intraAxonalVolumeImage->SetSpacing( m_WorkingSpacing ); intraAxonalVolumeImage->SetOrigin( m_WorkingOrigin ); intraAxonalVolumeImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); intraAxonalVolumeImage->SetLargestPossibleRegion( m_WorkingImageRegion ); intraAxonalVolumeImage->SetBufferedRegion( m_WorkingImageRegion ); intraAxonalVolumeImage->SetRequestedRegion( m_WorkingImageRegion ); intraAxonalVolumeImage->Allocate(); intraAxonalVolumeImage->FillBuffer(0); // determine intra-axonal volume fraction using the tract density itk::TractDensityImageFilter< ItkDoubleImgType >::Pointer tdiFilter = itk::TractDensityImageFilter< ItkDoubleImgType >::New(); tdiFilter->SetFiberBundle(m_FiberBundleTransformed); tdiFilter->SetBinaryOutput(false); tdiFilter->SetOutputAbsoluteValues(false); tdiFilter->SetInputImage(intraAxonalVolumeImage); tdiFilter->SetUseImageGeometry(true); tdiFilter->Update(); intraAxonalVolumeImage = tdiFilter->GetOutput(); m_StatusText += "0% 10 20 30 40 50 60 70 80 90 100%\n"; m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*"; boost::progress_display disp(m_TransformedMaskImage->GetLargestPossibleRegion().GetNumberOfPixels()*m_Parameters.m_SignalGen.GetNumVolumes()); for (unsigned int g=0; gSetSeed(signalModelSeed); for (int i=0; iSetSeed(signalModelSeed); if (m_Parameters.m_SignalGen.m_DoAddMotion && g>0) // if fibers have moved we need a new TDI and new directions { fOdfFilter->SetFiberBundle(m_FiberBundleTransformed); fOdfFilter->SetMaskImage(m_TransformedMaskImage); fOdfFilter->Update(); directionImageContainer = fOdfFilter->GetDirectionImageContainer(); tdiFilter->SetFiberBundle(m_FiberBundleTransformed); tdiFilter->Update(); intraAxonalVolumeImage = tdiFilter->GetOutput(); } ImageRegionIterator< ItkUcharImgType > it(m_TransformedMaskImage, m_TransformedMaskImage->GetLargestPossibleRegion()); while(!it.IsAtEnd()) { ++disp; unsigned long newTick = 50*disp.count()/disp.expected_count(); for (unsigned int tick = 0; tick<(newTick-lastTick); tick++) m_StatusText += "*"; lastTick = newTick; if (this->GetAbortGenerateData()) { m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n"; return; } if (it.Get()>0) { // generate fiber signal for (int c=0; cGetPixel(it.GetIndex()); for (unsigned int i=0; iSize(); i++) { itk::Vector< double, 3> dir; dir.CastFrom(directionImageContainer->GetElement(i)->GetPixel(it.GetIndex())); double norm = dir.GetNorm(); if (norm>0.0001) { m_Parameters.m_FiberModelList.at(c)->SetFiberDirection(dir); pix[g] += m_Parameters.m_FiberModelList.at(c)->SimulateMeasurement(g)*norm; count++; } } if (count>0) pix[g] /= count; pix[g] *= intraAxonalVolumeImage->GetPixel(it.GetIndex())*m_VoxelVolume; m_CompartmentImages.at(c)->SetPixel(it.GetIndex(), pix); } // simulate other compartments SimulateExtraAxonalSignal(it.GetIndex(), intraAxonalVolumeImage->GetPixel(it.GetIndex())*m_VoxelVolume, g); } ++it; } SimulateMotion(g); } itk::ImageFileWriter< ItkUcharImgType >::Pointer wr = itk::ImageFileWriter< ItkUcharImgType >::New(); wr->SetInput(fOdfFilter->GetNumDirectionsImage()); wr->SetFileName(mitk::IOUtil::GetTempPath()+"/NumDirections_MainFiberDirections.nrrd"); wr->Update(); break; } case (SignalGenerationParameters::RANDOM_DIRECTIONS): { ItkUcharImgType::Pointer numDirectionsImage = ItkUcharImgType::New(); numDirectionsImage->SetSpacing( m_WorkingSpacing ); numDirectionsImage->SetOrigin( m_WorkingOrigin ); numDirectionsImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); numDirectionsImage->SetLargestPossibleRegion( m_WorkingImageRegion ); numDirectionsImage->SetBufferedRegion( m_WorkingImageRegion ); numDirectionsImage->SetRequestedRegion( m_WorkingImageRegion ); numDirectionsImage->Allocate(); numDirectionsImage->FillBuffer(0); double sepAngle = cos(m_Parameters.m_SignalGen.m_FiberSeparationThreshold*M_PI/180.0); m_StatusText += "0% 10 20 30 40 50 60 70 80 90 100%\n"; m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*"; boost::progress_display disp(m_TransformedMaskImage->GetLargestPossibleRegion().GetNumberOfPixels()); ImageRegionIterator it(m_TransformedMaskImage, m_TransformedMaskImage->GetLargestPossibleRegion()); while(!it.IsAtEnd()) { ++disp; unsigned long newTick = 50*disp.count()/disp.expected_count(); for (unsigned int tick = 0; tick<(newTick-lastTick); tick++) m_StatusText += "*"; lastTick = newTick; if (this->GetAbortGenerateData()) { m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n"; return; } if (it.Get()>0) { int numFibs = m_RandGen->GetIntegerVariate(2)+1; DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(it.GetIndex()); double volume = m_RandGen->GetVariateWithClosedRange(0.3); // double sum = 0; std::vector< double > fractions; for (int i=0; iGetVariateWithClosedRange(0.5)); // sum += fractions.at(i); } // for (int i=0; i > directions; for (int i=0; iGetVariateWithClosedRange(2)-1.0; fib[1] = m_RandGen->GetVariateWithClosedRange(2)-1.0; fib[2] = m_RandGen->GetVariateWithClosedRange(2)-1.0; fib.Normalize(); double min = 0; for (unsigned int d=0; dmin) min = angle; } if (minSetFiberDirection(fib); pix += m_Parameters.m_FiberModelList.at(0)->SimulateMeasurement()*fractions[i]; directions.push_back(fib); } else i--; } pix *= (1-volume); m_CompartmentImages.at(0)->SetPixel(it.GetIndex(), pix); // CSF/GM { pix += volume*m_Parameters.m_NonFiberModelList.at(0)->SimulateMeasurement(); } numDirectionsImage->SetPixel(it.GetIndex(), numFibs); } ++it; } itk::ImageFileWriter< ItkUcharImgType >::Pointer wr = itk::ImageFileWriter< ItkUcharImgType >::New(); wr->SetInput(numDirectionsImage); wr->SetFileName(mitk::IOUtil::GetTempPath()+"/NumDirections_RandomDirections.nrrd"); wr->Update(); } } if (m_Logfile.is_open()) { m_Logfile << "DONE"; m_Logfile.close(); } m_StatusText += "\n\n"; if (this->GetAbortGenerateData()) { m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n"; return; } DoubleDwiType::Pointer doubleOutImage; double signalScale = m_Parameters.m_SignalGen.m_SignalScale; if ( m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition ) // do k-space stuff { m_StatusText += this->GetTime()+" > Simulating k-space acquisition using "+boost::lexical_cast(m_Parameters.m_SignalGen.m_NumberOfCoils)+" coil(s)\n"; MITK_INFO << "Simulating k-space acquisition using " << m_Parameters.m_SignalGen.m_NumberOfCoils << " coil(s)."; if (m_Parameters.m_SignalGen.m_DoSimulateRelaxation) m_StatusText += "Simulating signal relaxation\n"; if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull()) m_StatusText += "Simulating distortions\n"; if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging) m_StatusText += "Simulating ringing artifacts\n"; if (m_Parameters.m_SignalGen.m_EddyStrength>0) m_StatusText += "Simulating eddy currents\n"; if (m_Parameters.m_SignalGen.m_Spikes>0) m_StatusText += "Simulating spikes\n"; if (m_Parameters.m_SignalGen.m_CroppingFactor<1.0) m_StatusText += "Simulating aliasing artifacts\n"; if (m_Parameters.m_SignalGen.m_KspaceLineOffset>0) m_StatusText += "Simulating ghosts\n"; doubleOutImage = SimulateKspaceAcquisition(m_CompartmentImages); signalScale = 1; // already scaled in DoKspaceStuff } else // don't do k-space stuff, just sum compartments { m_StatusText += this->GetTime()+" > Summing compartments\n"; MITK_INFO << "Summing compartments"; doubleOutImage = m_CompartmentImages.at(0); for (unsigned int i=1; i::Pointer adder = itk::AddImageFilter< DoubleDwiType, DoubleDwiType, DoubleDwiType>::New(); adder->SetInput1(doubleOutImage); adder->SetInput2(m_CompartmentImages.at(i)); adder->Update(); doubleOutImage = adder->GetOutput(); } } if (this->GetAbortGenerateData()) { m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n"; return; } m_StatusText += this->GetTime()+" > Finalizing image\n"; MITK_INFO << "Finalizing image"; if (signalScale>1) m_StatusText += " Scaling signal\n"; if (m_Parameters.m_NoiseModel!=NULL) m_StatusText += " Adding noise\n"; unsigned int window = 0; unsigned int min = itk::NumericTraits::max(); ImageRegionIterator it4 (m_OutputImage, m_OutputImage->GetLargestPossibleRegion()); DoubleDwiType::PixelType signal; signal.SetSize(m_Parameters.m_SignalGen.GetNumVolumes()); boost::progress_display disp2(m_OutputImage->GetLargestPossibleRegion().GetNumberOfPixels()); m_StatusText += "0% 10 20 30 40 50 60 70 80 90 100%\n"; m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*"; lastTick = 0; while(!it4.IsAtEnd()) { if (this->GetAbortGenerateData()) { m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n"; return; } ++disp2; unsigned long newTick = 50*disp2.count()/disp2.expected_count(); for (unsigned long tick = 0; tick<(newTick-lastTick); tick++) m_StatusText += "*"; lastTick = newTick; typename OutputImageType::IndexType index = it4.GetIndex(); signal = doubleOutImage->GetPixel(index)*signalScale; if (m_Parameters.m_NoiseModel!=NULL) m_Parameters.m_NoiseModel->AddNoise(signal); for (unsigned int i=0; i0) signal[i] = floor(signal[i]+0.5); else signal[i] = ceil(signal[i]-0.5); if ( (!m_Parameters.m_SignalGen.IsBaselineIndex(i) || signal.Size()==1) && signal[i]>window) window = signal[i]; if ( (!m_Parameters.m_SignalGen.IsBaselineIndex(i) || signal.Size()==1) && signal[i]SetNthOutput(0, m_OutputImage); m_StatusText += "\n\n"; m_StatusText += "Finished simulation\n"; m_StatusText += "Simulation time: "+GetTime(); m_TimeProbe.Stop(); } template< class PixelType > void TractsToDWIImageFilter< PixelType >::SimulateMotion(int g) { if (m_Parameters.m_SignalGen.m_DoAddMotion && gGetDeepCopy(); m_Rotation[0] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[0]*2)-m_Parameters.m_SignalGen.m_Rotation[0]; m_Rotation[1] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[1]*2)-m_Parameters.m_SignalGen.m_Rotation[1]; m_Rotation[2] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[2]*2)-m_Parameters.m_SignalGen.m_Rotation[2]; m_Translation[0] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[0]*2)-m_Parameters.m_SignalGen.m_Translation[0]; m_Translation[1] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[1]*2)-m_Parameters.m_SignalGen.m_Translation[1]; m_Translation[2] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[2]*2)-m_Parameters.m_SignalGen.m_Translation[2]; } // rotate mask image if (m_MaskImageSet) { ImageRegionIterator maskIt(m_UpsampledMaskImage, m_UpsampledMaskImage->GetLargestPossibleRegion()); m_TransformedMaskImage->FillBuffer(0); while(!maskIt.IsAtEnd()) { if (maskIt.Get()<=0) { ++maskIt; continue; } DoubleDwiType::IndexType index = maskIt.GetIndex(); itk::Point point; m_UpsampledMaskImage->TransformIndexToPhysicalPoint(index, point); if (m_Parameters.m_SignalGen.m_DoRandomizeMotion) point = m_FiberBundleWorkingCopy->TransformPoint(point.GetVnlVector(), m_Rotation[0],m_Rotation[1],m_Rotation[2],m_Translation[0],m_Translation[1],m_Translation[2]); else point = m_FiberBundleWorkingCopy->TransformPoint(point.GetVnlVector(), m_Rotation[0]*(g+1),m_Rotation[1]*(g+1),m_Rotation[2]*(g+1),m_Translation[0]*(g+1),m_Translation[1]*(g+1),m_Translation[2]*(g+1)); m_TransformedMaskImage->TransformPhysicalPointToIndex(point, index); if (m_TransformedMaskImage->GetLargestPossibleRegion().IsInside(index)) m_TransformedMaskImage->SetPixel(index,100); ++maskIt; } } if (m_Parameters.m_SignalGen.m_DoRandomizeMotion) { m_Rotations.push_back(m_Rotation); m_Translations.push_back(m_Translation); } else { m_Rotations.push_back(m_Rotation*(g+1)); m_Translations.push_back(m_Translation*(g+1)); } // rotate fibers if (m_Logfile.is_open()) { m_Logfile << g+1 << " rotation: " << m_Rotation[0] << "," << m_Rotation[1] << "," << m_Rotation[2] << ";"; m_Logfile << " translation: " << m_Translation[0] << "," << m_Translation[1] << "," << m_Translation[2] << "\n"; } m_FiberBundleTransformed->TransformFibers(m_Rotation[0],m_Rotation[1],m_Rotation[2],m_Translation[0],m_Translation[1],m_Translation[2]); } } template< class PixelType > void TractsToDWIImageFilter< PixelType >::SimulateExtraAxonalSignal(ItkUcharImgType::IndexType index, double intraAxonalVolume, int g) { int numFiberCompartments = m_Parameters.m_FiberModelList.size(); int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size(); if (intraAxonalVolume>0.0001 && m_Parameters.m_SignalGen.m_DoDisablePartialVolume) // only fiber in voxel { DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(index); if (g>=0) pix[g] *= m_VoxelVolume/intraAxonalVolume; else pix *= m_VoxelVolume/intraAxonalVolume; m_CompartmentImages.at(0)->SetPixel(index, pix); if (g==0) m_VolumeFractions.at(0)->SetPixel(index, 1); for (int i=1; iGetPixel(index); if (g>=0) pix[g] = 0.0; else pix.Fill(0.0); m_CompartmentImages.at(i)->SetPixel(index, pix); } } else { if (g==0) m_VolumeFractions.at(0)->SetPixel(index, intraAxonalVolume/m_VoxelVolume); // get non-transformed point (remove headmotion tranformation) // this point can then be transformed to each of the original images, regardless of their geometry itk::Point point; m_TransformedMaskImage->TransformIndexToPhysicalPoint(index, point); if (m_Parameters.m_SignalGen.m_DoAddMotion) { if (m_Parameters.m_SignalGen.m_DoRandomizeMotion && g>0) point = m_FiberBundleWorkingCopy->TransformPoint(point.GetVnlVector(), -m_Rotation[0],-m_Rotation[1],-m_Rotation[2],-m_Translation[0],-m_Translation[1],-m_Translation[2]); else if (g>=0) point = m_FiberBundleWorkingCopy->TransformPoint(point.GetVnlVector(), -m_Rotation[0]*g,-m_Rotation[1]*g,-m_Rotation[2]*g,-m_Translation[0]*g,-m_Translation[1]*g,-m_Translation[2]*g); } if (m_Parameters.m_SignalGen.m_DoDisablePartialVolume) { int maxVolumeIndex = 0; double maxWeight = 0; for (int i=0; i1) { // ItkUcharImgType::IndexType maskIndex; // m_UpsampledMaskImage->TransformPhysicalPointToIndex(point, maskIndex); // if (!m_UpsampledMaskImage->GetLargestPossibleRegion().IsInside(maskIndex) || m_UpsampledMaskImage->GetPixel(maskIndex)==0) // continue; DoubleDwiType::IndexType newIndex; m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->TransformPhysicalPointToIndex(point, newIndex); if (!m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetLargestPossibleRegion().IsInside(newIndex)) continue; weight = m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetPixel(newIndex); } if (weight>maxWeight) { maxWeight = weight; maxVolumeIndex = i; } } DoubleDwiType::Pointer doubleDwi = m_CompartmentImages.at(maxVolumeIndex+numFiberCompartments); DoubleDwiType::PixelType pix = doubleDwi->GetPixel(index); if (g>=0) pix[g] += m_Parameters.m_NonFiberModelList[maxVolumeIndex]->SimulateMeasurement(g); else pix += m_Parameters.m_NonFiberModelList[maxVolumeIndex]->SimulateMeasurement(); doubleDwi->SetPixel(index, pix); if (g==0) m_VolumeFractions.at(maxVolumeIndex+numFiberCompartments)->SetPixel(index, 1); } else { double extraAxonalVolume = m_VoxelVolume-intraAxonalVolume; // non-fiber volume if (extraAxonalVolume<0) { MITK_ERROR << "Coorupted intra-axonal signal voxel detected. Fiber volume larger voxel volume!"; extraAxonalVolume = 0; } double interAxonalVolume = 0; if (numFiberCompartments>1) interAxonalVolume = extraAxonalVolume * intraAxonalVolume/m_VoxelVolume; // inter-axonal fraction of non fiber compartment double other = extraAxonalVolume - interAxonalVolume; // rest of compartment if (other<0) { MITK_ERROR << "Coorupted signal voxel detected. Fiber volume larger voxel volume!"; other = 0; } // adjust non-fiber and intra-axonal signal for (int i=1; iGetPixel(index); if (intraAxonalVolume>0) // remove scaling by intra-axonal volume from inter-axonal compartment { if (g>=0) pix[g] /= intraAxonalVolume; else pix /= intraAxonalVolume; } if (m_Parameters.m_FiberModelList[i]->GetVolumeFractionImage()!=nullptr) { // ItkUcharImgType::IndexType maskIndex; // m_UpsampledMaskImage->TransformPhysicalPointToIndex(point, maskIndex); // if (!m_UpsampledMaskImage->GetLargestPossibleRegion().IsInside(maskIndex) || m_UpsampledMaskImage->GetPixel(maskIndex)==0) // continue; DoubleDwiType::IndexType newIndex; m_Parameters.m_FiberModelList[i]->GetVolumeFractionImage()->TransformPhysicalPointToIndex(point, newIndex); if (!m_Parameters.m_FiberModelList[i]->GetVolumeFractionImage()->GetLargestPossibleRegion().IsInside(newIndex)) continue; weight = m_Parameters.m_FiberModelList[i]->GetVolumeFractionImage()->GetPixel(newIndex)*m_VoxelVolume; } if (g>=0) pix[g] *= weight; else pix *= weight; m_CompartmentImages.at(i)->SetPixel(index, pix); if (g==0) m_VolumeFractions.at(i)->SetPixel(index, weight/m_VoxelVolume); } for (int i=0; iGetPixel(index); if (m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()!=nullptr) { // ItkUcharImgType::IndexType maskIndex; // m_UpsampledMaskImage->TransformPhysicalPointToIndex(point, maskIndex); // if (!m_UpsampledMaskImage->GetLargestPossibleRegion().IsInside(maskIndex) || m_UpsampledMaskImage->GetPixel(maskIndex)==0) // continue; DoubleDwiType::IndexType newIndex; m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->TransformPhysicalPointToIndex(point, newIndex); if (!m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetLargestPossibleRegion().IsInside(newIndex)) continue; weight = m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetPixel(newIndex)*m_VoxelVolume; if (m_UseRelativeNonFiberVolumeFractions) weight *= other/m_VoxelVolume; } if (g>=0) pix[g] += m_Parameters.m_NonFiberModelList[i]->SimulateMeasurement(g)*weight; else pix += m_Parameters.m_NonFiberModelList[i]->SimulateMeasurement()*weight; m_CompartmentImages.at(i+numFiberCompartments)->SetPixel(index, pix); if (g==0) m_VolumeFractions.at(i+numFiberCompartments)->SetPixel(index, weight/m_VoxelVolume); } } } } template< class PixelType > itk::Point TractsToDWIImageFilter< PixelType >::GetItkPoint(double point[3]) { itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; return itkPoint; } template< class PixelType > itk::Vector TractsToDWIImageFilter< PixelType >::GetItkVector(double point[3]) { itk::Vector itkVector; itkVector[0] = point[0]; itkVector[1] = point[1]; itkVector[2] = point[2]; return itkVector; } template< class PixelType > vnl_vector_fixed TractsToDWIImageFilter< PixelType >::GetVnlVector(double point[3]) { vnl_vector_fixed vnlVector; vnlVector[0] = point[0]; vnlVector[1] = point[1]; vnlVector[2] = point[2]; return vnlVector; } template< class PixelType > vnl_vector_fixed TractsToDWIImageFilter< PixelType >::GetVnlVector(Vector& vector) { vnl_vector_fixed vnlVector; vnlVector[0] = vector[0]; vnlVector[1] = vector[1]; vnlVector[2] = vector[2]; return vnlVector; } template< class PixelType > double TractsToDWIImageFilter< PixelType >::RoundToNearest(double num) { return (num > 0.0) ? floor(num + 0.5) : ceil(num - 0.5); } template< class PixelType > std::string TractsToDWIImageFilter< PixelType >::GetTime() { m_TimeProbe.Stop(); unsigned long total = RoundToNearest(m_TimeProbe.GetTotal()); unsigned long hours = total/3600; unsigned long minutes = (total%3600)/60; unsigned long seconds = total%60; std::string out = ""; out.append(boost::lexical_cast(hours)); out.append(":"); out.append(boost::lexical_cast(minutes)); out.append(":"); out.append(boost::lexical_cast(seconds)); m_TimeProbe.Start(); return out; } } diff --git a/Modules/DiffusionImaging/MiniApps/PeakExtraction.cpp b/Modules/DiffusionImaging/MiniApps/PeakExtraction.cpp index 17d8b4764f..2a39c3dcb5 100755 --- a/Modules/DiffusionImaging/MiniApps/PeakExtraction.cpp +++ b/Modules/DiffusionImaging/MiniApps/PeakExtraction.cpp @@ -1,374 +1,380 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include "mitkCommandLineParser.h" #include #include #include #include #include template int StartPeakExtraction(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); parser.addArgument("image", "i", mitkCommandLineParser::InputFile, "Input image", "sh coefficient image", us::Any(), false); parser.addArgument("outroot", "o", mitkCommandLineParser::OutputDirectory, "Output directory", "output root", us::Any(), false); parser.addArgument("mask", "m", mitkCommandLineParser::InputFile, "Mask", "mask image"); parser.addArgument("normalization", "n", mitkCommandLineParser::Int, "Normalization", "0=no norm, 1=max norm, 2=single vec norm", 1, true); parser.addArgument("numpeaks", "p", mitkCommandLineParser::Int, "Max. number of peaks", "maximum number of extracted peaks", 2, true); parser.addArgument("peakthres", "r", mitkCommandLineParser::Float, "Peak threshold", "peak threshold relative to largest peak", 0.4, true); parser.addArgument("abspeakthres", "a", mitkCommandLineParser::Float, "Absolute peak threshold", "absolute peak threshold weighted with local GFA value", 0.06, true); parser.addArgument("shConvention", "s", mitkCommandLineParser::String, "Use specified SH-basis", "use specified SH-basis (MITK, FSL, MRtrix)", string("MITK"), true); parser.addArgument("noFlip", "f", mitkCommandLineParser::Bool, "No flip", "do not flip input image to match MITK coordinate convention"); + parser.addArgument("clusterThres", "c", mitkCommandLineParser::Float, "Clustering threshold", "directions closer together than the specified angular threshold will be clustered (in rad)", 0.9); parser.addArgument("flipX", "fx", mitkCommandLineParser::Bool, "Flip X", "Flip peaks in x direction"); parser.addArgument("flipY", "fy", mitkCommandLineParser::Bool, "Flip Y", "Flip peaks in y direction"); parser.addArgument("flipZ", "fz", mitkCommandLineParser::Bool, "Flip Z", "Flip peaks in z direction"); parser.setCategory("Preprocessing Tools"); parser.setTitle("Peak Extraction"); parser.setDescription(""); parser.setContributor("MBI"); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments string imageName = us::any_cast(parsedArgs["image"]); string outRoot = us::any_cast(parsedArgs["outroot"]); // optional arguments string maskImageName(""); if (parsedArgs.count("mask")) maskImageName = us::any_cast(parsedArgs["mask"]); int normalization = 1; if (parsedArgs.count("normalization")) normalization = us::any_cast(parsedArgs["normalization"]); int numPeaks = 2; if (parsedArgs.count("numpeaks")) numPeaks = us::any_cast(parsedArgs["numpeaks"]); float peakThres = 0.4; if (parsedArgs.count("peakthres")) peakThres = us::any_cast(parsedArgs["peakthres"]); float absPeakThres = 0.06; if (parsedArgs.count("abspeakthres")) absPeakThres = us::any_cast(parsedArgs["abspeakthres"]); + float clusterThres = 0.9; + if (parsedArgs.count("clusterThres")) + clusterThres = us::any_cast(parsedArgs["clusterThres"]); + bool noFlip = false; if (parsedArgs.count("noFlip")) noFlip = us::any_cast(parsedArgs["noFlip"]); bool flipX = false; if (parsedArgs.count("flipX")) flipX = us::any_cast(parsedArgs["flipX"]); bool flipY = false; if (parsedArgs.count("flipY")) flipY = us::any_cast(parsedArgs["flipY"]); bool flipZ = false; if (parsedArgs.count("flipZ")) flipZ = us::any_cast(parsedArgs["flipZ"]); std::cout << "image: " << imageName; std::cout << "outroot: " << outRoot; if (!maskImageName.empty()) std::cout << "mask: " << maskImageName; else std::cout << "no mask image selected"; std::cout << "numpeaks: " << numPeaks; std::cout << "peakthres: " << peakThres; std::cout << "abspeakthres: " << absPeakThres; std::cout << "shOrder: " << shOrder; try { mitk::Image::Pointer image = mitk::IOUtil::LoadImage(imageName); mitk::Image::Pointer mask = mitk::IOUtil::LoadImage(maskImageName); typedef itk::Image ItkUcharImgType; typedef itk::FiniteDiffOdfMaximaExtractionFilter< float, shOrder, 20242 > MaximaExtractionFilterType; typename MaximaExtractionFilterType::Pointer filter = MaximaExtractionFilterType::New(); int toolkitConvention = 0; if (parsedArgs.count("shConvention")) { string convention = us::any_cast(parsedArgs["shConvention"]).c_str(); if ( boost::algorithm::equals(convention, "FSL") ) { toolkitConvention = 1; std::cout << "Using FSL SH-basis"; } else if ( boost::algorithm::equals(convention, "MRtrix") ) { toolkitConvention = 2; std::cout << "Using MRtrix SH-basis"; } else std::cout << "Using MITK SH-basis"; } else std::cout << "Using MITK SH-basis"; ItkUcharImgType::Pointer itkMaskImage = NULL; if (mask.IsNotNull()) { try{ itkMaskImage = ItkUcharImgType::New(); mitk::CastToItkImage(mask, itkMaskImage); filter->SetMaskImage(itkMaskImage); } catch(...) { } } if (toolkitConvention>0) { std::cout << "Converting coefficient image to MITK format"; typedef itk::ShCoefficientImageImporter< float, shOrder > ConverterType; typedef mitk::ImageToItk< itk::Image< float, 4 > > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(image); caster->Update(); itk::Image< float, 4 >::Pointer itkImage = caster->GetOutput(); typename ConverterType::Pointer converter = ConverterType::New(); if (noFlip) { converter->SetInputImage(itkImage); } else { std::cout << "Flipping image"; itk::FixedArray flipAxes; flipAxes[0] = true; flipAxes[1] = true; flipAxes[2] = false; flipAxes[3] = false; itk::FlipImageFilter< itk::Image< float, 4 > >::Pointer flipper = itk::FlipImageFilter< itk::Image< float, 4 > >::New(); flipper->SetInput(itkImage); flipper->SetFlipAxes(flipAxes); flipper->Update(); itk::Image< float, 4 >::Pointer flipped = flipper->GetOutput(); itk::Matrix< double,4,4 > m = itkImage->GetDirection(); m[0][0] *= -1; m[1][1] *= -1; flipped->SetDirection(m); itk::Point< float, 4 > o = itkImage->GetOrigin(); o[0] -= (flipped->GetLargestPossibleRegion().GetSize(0)-1); o[1] -= (flipped->GetLargestPossibleRegion().GetSize(1)-1); flipped->SetOrigin(o); converter->SetInputImage(flipped); } std::cout << "Starting conversion"; switch (toolkitConvention) { case 1: converter->SetToolkit(ConverterType::FSL); filter->SetToolkit(MaximaExtractionFilterType::FSL); break; case 2: converter->SetToolkit(ConverterType::MRTRIX); filter->SetToolkit(MaximaExtractionFilterType::MRTRIX); break; default: converter->SetToolkit(ConverterType::FSL); filter->SetToolkit(MaximaExtractionFilterType::FSL); break; } converter->GenerateData(); filter->SetInput(converter->GetCoefficientImage()); } else { try{ typedef mitk::ImageToItk< typename MaximaExtractionFilterType::CoefficientImageType > CasterType; typename CasterType::Pointer caster = CasterType::New(); caster->SetInput(image); caster->Update(); filter->SetInput(caster->GetOutput()); } catch(...) { std::cout << "wrong image type"; return EXIT_FAILURE; } } filter->SetMaxNumPeaks(numPeaks); filter->SetPeakThreshold(peakThres); filter->SetAbsolutePeakThreshold(absPeakThres); filter->SetAngularThreshold(1); + filter->SetClusteringThreshold(clusterThres); filter->SetFlipX(flipX); filter->SetFlipY(flipY); filter->SetFlipZ(flipZ); switch (normalization) { case 0: filter->SetNormalizationMethod(MaximaExtractionFilterType::NO_NORM); break; case 1: filter->SetNormalizationMethod(MaximaExtractionFilterType::MAX_VEC_NORM); break; case 2: filter->SetNormalizationMethod(MaximaExtractionFilterType::SINGLE_VEC_NORM); break; } std::cout << "Starting extraction"; filter->Update(); // write direction images { typedef typename MaximaExtractionFilterType::ItkDirectionImageContainer ItkDirectionImageContainer; typename ItkDirectionImageContainer::Pointer container = filter->GetDirectionImageContainer(); for (unsigned int i=0; iSize(); i++) { typename MaximaExtractionFilterType::ItkDirectionImage::Pointer itkImg = container->GetElement(i); if (itkMaskImage.IsNotNull()) { itkImg->SetDirection(itkMaskImage->GetDirection()); itkImg->SetOrigin(itkMaskImage->GetOrigin()); } string outfilename = outRoot; outfilename.append("_DIRECTION_"); outfilename.append(boost::lexical_cast(i)); outfilename.append(".nrrd"); typedef itk::ImageFileWriter< typename MaximaExtractionFilterType::ItkDirectionImage > WriterType; typename WriterType::Pointer writer = WriterType::New(); writer->SetFileName(outfilename); writer->SetInput(itkImg); writer->Update(); } } // write num directions image { ItkUcharImgType::Pointer numDirImage = filter->GetNumDirectionsImage(); if (itkMaskImage.IsNotNull()) { numDirImage->SetDirection(itkMaskImage->GetDirection()); numDirImage->SetOrigin(itkMaskImage->GetOrigin()); } string outfilename = outRoot.c_str(); outfilename.append("_NUM_DIRECTIONS.nrrd"); typedef itk::ImageFileWriter< ItkUcharImgType > WriterType; WriterType::Pointer writer = WriterType::New(); writer->SetFileName(outfilename); writer->SetInput(numDirImage); writer->Update(); } // write vector field { mitk::FiberBundle::Pointer directions = filter->GetOutputFiberBundle(); string outfilename = outRoot.c_str(); outfilename.append("_VECTOR_FIELD.fib"); mitk::IOUtil::Save(directions.GetPointer(),outfilename.c_str()); } } catch (itk::ExceptionObject e) { std::cout << e; return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); parser.addArgument("image", "i", mitkCommandLineParser::InputFile, "Input image", "sh coefficient image", us::Any(), false); parser.addArgument("shOrder", "sh", mitkCommandLineParser::Int, "Spherical harmonics order", "spherical harmonics order"); parser.addArgument("outroot", "o", mitkCommandLineParser::OutputDirectory, "Output directory", "output root", us::Any(), false); parser.addArgument("mask", "m", mitkCommandLineParser::InputFile, "Mask", "mask image"); parser.addArgument("normalization", "n", mitkCommandLineParser::Int, "Normalization", "0=no norm, 1=max norm, 2=single vec norm", 1, true); parser.addArgument("numpeaks", "p", mitkCommandLineParser::Int, "Max. number of peaks", "maximum number of extracted peaks", 2, true); parser.addArgument("peakthres", "r", mitkCommandLineParser::Float, "Peak threshold", "peak threshold relative to largest peak", 0.4, true); parser.addArgument("abspeakthres", "a", mitkCommandLineParser::Float, "Absolute peak threshold", "absolute peak threshold weighted with local GFA value", 0.06, true); parser.addArgument("shConvention", "s", mitkCommandLineParser::String, "Use specified SH-basis", "use specified SH-basis (MITK, FSL, MRtrix)", string("MITK"), true); parser.addArgument("noFlip", "f", mitkCommandLineParser::Bool, "No flip", "do not flip input image to match MITK coordinate convention"); parser.setCategory("Preprocessing Tools"); parser.setTitle("Peak Extraction"); parser.setDescription(""); parser.setContributor("MBI"); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; int shOrder = -1; if (parsedArgs.count("shOrder")) shOrder = us::any_cast(parsedArgs["shOrder"]); switch (shOrder) { case 4: return StartPeakExtraction<4>(argc, argv); case 6: return StartPeakExtraction<6>(argc, argv); case 8: return StartPeakExtraction<8>(argc, argv); case 10: return StartPeakExtraction<10>(argc, argv); case 12: return StartPeakExtraction<12>(argc, argv); } return EXIT_FAILURE; }