diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkAddArtifactsToDwiImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkAddArtifactsToDwiImageFilter.cpp index a8c4330d70..4e9c2430d8 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkAddArtifactsToDwiImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkAddArtifactsToDwiImageFilter.cpp @@ -1,336 +1,336 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkAddArtifactsToDwiImageFilter_txx #define __itkAddArtifactsToDwiImageFilter_txx #include #include #include #include "itkAddArtifactsToDwiImageFilter.h" #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include namespace itk { template< class TPixelType > AddArtifactsToDwiImageFilter< TPixelType > ::AddArtifactsToDwiImageFilter() : m_UseConstantRandSeed(false) { this->SetNumberOfRequiredInputs( 1 ); m_RandGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); m_RandGen->SetSeed(); } template< class TPixelType > void AddArtifactsToDwiImageFilter< TPixelType > ::GenerateData() { if (m_UseConstantRandSeed) // always generate the same random numbers? m_RandGen->SetSeed(0); else m_RandGen->SetSeed(); m_StartTime = clock(); m_StatusText = "Starting simulation\n"; typename InputImageType::Pointer inputImage = static_cast< InputImageType* >( this->ProcessObject::GetInput(0) ); itk::ImageRegion<3> inputRegion = inputImage->GetLargestPossibleRegion(); typename itk::ImageDuplicator::Pointer duplicator = itk::ImageDuplicator::New(); duplicator->SetInputImage( inputImage ); duplicator->Update(); typename InputImageType::Pointer outputImage = duplicator->GetOutput(); // is input slize size even? int xMax=inputRegion.GetSize(0); int yMax=inputRegion.GetSize(1); if ( xMax%2 == 1 ) xMax += 1; if ( yMax%2 == 1 ) yMax += 1; // create slice object typename SliceType::Pointer slice = SliceType::New(); ImageRegion<2> sliceRegion; sliceRegion.SetSize(0, xMax); sliceRegion.SetSize(1, yMax); slice->SetLargestPossibleRegion( sliceRegion ); slice->SetBufferedRegion( sliceRegion ); slice->SetRequestedRegion( sliceRegion ); slice->Allocate(); slice->FillBuffer(0.0); ImageRegion<2> upsampledSliceRegion; - if (m_Parameters.m_DoAddGibbsRinging) + if ( m_Parameters.m_SignalGen.m_DoAddGibbsRinging) { upsampledSliceRegion.SetSize(0, xMax*2); upsampledSliceRegion.SetSize(1, yMax*2); } // frequency map slice typename SliceType::Pointer fMapSlice = NULL; - if (m_Parameters.m_FrequencyMap.IsNotNull()) + if ( m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull()) { fMapSlice = SliceType::New(); fMapSlice->SetLargestPossibleRegion( sliceRegion ); fMapSlice->SetBufferedRegion( sliceRegion ); fMapSlice->SetRequestedRegion( sliceRegion ); fMapSlice->Allocate(); fMapSlice->FillBuffer(0.0); } - m_Parameters.m_SignalScale = 1; - m_Parameters.m_DoSimulateRelaxation = false; + m_Parameters.m_SignalGen.m_SignalScale = 1; + m_Parameters.m_SignalGen.m_DoSimulateRelaxation = false; - if (m_Parameters.m_Spikes>0 || m_Parameters.m_FrequencyMap.IsNotNull() || m_Parameters.m_KspaceLineOffset>0.0 || m_Parameters.m_DoAddGibbsRinging || m_Parameters.m_EddyStrength>0 || m_Parameters.m_CroppingFactor<1.0) + if ( m_Parameters.m_SignalGen.m_Spikes>0 || m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull() || m_Parameters.m_SignalGen.m_KspaceLineOffset>0.0 || m_Parameters.m_SignalGen.m_DoAddGibbsRinging || m_Parameters.m_SignalGen.m_EddyStrength>0 || m_Parameters.m_SignalGen.m_CroppingFactor<1.0) { - ImageRegion<3> croppedRegion = inputRegion; croppedRegion.SetSize(1, croppedRegion.GetSize(1)*m_Parameters.m_CroppingFactor); + ImageRegion<3> croppedRegion = inputRegion; croppedRegion.SetSize(1, croppedRegion.GetSize(1)* m_Parameters.m_SignalGen.m_CroppingFactor); itk::Point shiftedOrigin = inputImage->GetOrigin(); shiftedOrigin[1] += (inputRegion.GetSize(1)-croppedRegion.GetSize(1))*inputImage->GetSpacing()[1]/2; outputImage = InputImageType::New(); outputImage->SetSpacing( inputImage->GetSpacing() ); outputImage->SetOrigin( shiftedOrigin ); outputImage->SetDirection( inputImage->GetDirection() ); outputImage->SetLargestPossibleRegion( croppedRegion ); outputImage->SetBufferedRegion( croppedRegion ); outputImage->SetRequestedRegion( croppedRegion ); outputImage->SetVectorLength( inputImage->GetVectorLength() ); outputImage->Allocate(); typename InputImageType::PixelType temp; temp.SetSize(inputImage->GetVectorLength()); temp.Fill(0.0); outputImage->FillBuffer(temp); int tempY=croppedRegion.GetSize(1); tempY += tempY%2; croppedRegion.SetSize(1, tempY); m_StatusText += this->GetTime()+" > Adjusting complex signal\n"; - if (m_Parameters.m_FrequencyMap.IsNotNull()) + if ( m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull()) m_StatusText += "Simulating distortions\n"; - if (m_Parameters.m_DoAddGibbsRinging) + if ( m_Parameters.m_SignalGen.m_DoAddGibbsRinging) m_StatusText += "Simulating ringing artifacts\n"; - if (m_Parameters.m_EddyStrength>0) + if ( m_Parameters.m_SignalGen.m_EddyStrength>0) m_StatusText += "Simulating eddy currents\n"; - if (m_Parameters.m_Spikes>0) + if ( m_Parameters.m_SignalGen.m_Spikes>0) m_StatusText += "Simulating spikes\n"; - if (m_Parameters.m_CroppingFactor<1.0) + if ( m_Parameters.m_SignalGen.m_CroppingFactor<1.0) m_StatusText += "Simulating aliasing artifacts\n"; - if (m_Parameters.m_KspaceLineOffset>0) + if ( m_Parameters.m_SignalGen.m_KspaceLineOffset>0) m_StatusText += "Simulating ghosts\n"; std::vector< unsigned int > spikeVolume; - for (unsigned int i=0; iGetIntegerVariate()%inputImage->GetVectorLength()); std::sort (spikeVolume.begin(), spikeVolume.end()); std::reverse (spikeVolume.begin(), spikeVolume.end()); - FiberfoxParameters doubleParam = m_Parameters.CopyParameters(); + FiberfoxParameters doubleParam = m_Parameters.CopyParameters(); m_StatusText += "0% 10 20 30 40 50 60 70 80 90 100%\n"; m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*"; unsigned long lastTick = 0; boost::progress_display disp(inputImage->GetVectorLength()*inputRegion.GetSize(2)); for (unsigned int g=0; gGetVectorLength(); g++) { std::vector< unsigned int > spikeSlice; while (!spikeVolume.empty() && spikeVolume.back()==g) { spikeSlice.push_back(m_RandGen->GetIntegerVariate()%inputImage->GetLargestPossibleRegion().GetSize(2)); spikeVolume.pop_back(); } std::sort (spikeSlice.begin(), spikeSlice.end()); std::reverse (spikeSlice.begin(), spikeSlice.end()); for (unsigned int z=0; zGetAbortGenerateData()) { m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n"; return; } std::vector< SliceType::Pointer > compartmentSlices; // extract slice from channel g for (unsigned int y=0; yGetPixel(index3D)[g]; slice->SetPixel(index2D, pix2D); if (fMapSlice.IsNotNull()) - fMapSlice->SetPixel(index2D, m_Parameters.m_FrequencyMap->GetPixel(index3D)); + fMapSlice->SetPixel(index2D, m_Parameters.m_SignalGen.m_FrequencyMap->GetPixel(index3D)); } - if (m_Parameters.m_DoAddGibbsRinging) + if ( m_Parameters.m_SignalGen.m_DoAddGibbsRinging) { itk::ResampleImageFilter::Pointer resampler = itk::ResampleImageFilter::New(); resampler->SetInput(slice); resampler->SetOutputParametersFromImage(slice); resampler->SetSize(upsampledSliceRegion.GetSize()); resampler->SetOutputSpacing(slice->GetSpacing()/2); itk::NearestNeighborInterpolateImageFunction::Pointer nn_interpolator = itk::NearestNeighborInterpolateImageFunction::New(); resampler->SetInterpolator(nn_interpolator); resampler->Update(); typename SliceType::Pointer upslice = resampler->GetOutput(); compartmentSlices.push_back(upslice); if (fMapSlice.IsNotNull()) { itk::ResampleImageFilter::Pointer resampler = itk::ResampleImageFilter::New(); resampler->SetInput(fMapSlice); resampler->SetOutputParametersFromImage(fMapSlice); resampler->SetSize(upsampledSliceRegion.GetSize()); resampler->SetOutputSpacing(fMapSlice->GetSpacing()/2); itk::NearestNeighborInterpolateImageFunction::Pointer nn_interpolator = itk::NearestNeighborInterpolateImageFunction::New(); resampler->SetInterpolator(nn_interpolator); resampler->Update(); fMapSlice = resampler->GetOutput(); } } else compartmentSlices.push_back(slice); // fourier transform slice typename ComplexSliceType::Pointer fSlice; itk::Size<2> outSize; outSize.SetElement(0, xMax); outSize.SetElement(1, croppedRegion.GetSize()[1]); typename itk::KspaceImageFilter< SliceType::PixelType >::Pointer idft = itk::KspaceImageFilter< SliceType::PixelType >::New(); idft->SetUseConstantRandSeed(m_UseConstantRandSeed); idft->SetParameters(doubleParam); idft->SetCompartmentImages(compartmentSlices); idft->SetFrequencyMapSlice(fMapSlice); - idft->SetDiffusionGradientDirection(m_Parameters.GetGradientDirection(g)); + idft->SetDiffusionGradientDirection( m_Parameters.GetGradientDirection(g)); idft->SetZ((double)z-(double)inputRegion.GetSize(2)/2.0); idft->SetOutSize(outSize); int numSpikes = 0; while (!spikeSlice.empty() && spikeSlice.back()==z) { numSpikes++; spikeSlice.pop_back(); } idft->SetSpikesPerSlice(numSpikes); idft->Update(); fSlice = idft->GetOutput(); // inverse fourier transform slice typename SliceType::Pointer newSlice; typename itk::DftImageFilter< SliceType::PixelType >::Pointer dft = itk::DftImageFilter< SliceType::PixelType >::New(); dft->SetInput(fSlice); dft->Update(); newSlice = dft->GetOutput(); // put slice back into channel g for (unsigned int y=0; yGetLargestPossibleRegion().GetSize(1); y++) for (unsigned int x=0; xGetLargestPossibleRegion().GetSize(0); x++) { typename InputImageType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z; typename InputImageType::PixelType pix3D = outputImage->GetPixel(index3D); typename SliceType::IndexType index2D; index2D[0]=x; index2D[1]=y; double signal = newSlice->GetPixel(index2D); if (signal>0) signal = floor(signal+0.5); else signal = ceil(signal-0.5); pix3D[g] = signal; outputImage->SetPixel(index3D, pix3D); } ++disp; unsigned long newTick = 50*disp.count()/disp.expected_count(); for (unsigned int tick = 0; tick<(newTick-lastTick); tick++) m_StatusText += "*"; lastTick = newTick; } } m_StatusText += "\n\n"; } - if (m_Parameters.m_NoiseModel!=NULL) + if ( m_Parameters.m_NoiseModel!=NULL) { m_StatusText += this->GetTime()+" > Adding noise\n"; m_StatusText += "0% 10 20 30 40 50 60 70 80 90 100%\n"; m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*"; unsigned long lastTick = 0; ImageRegionIterator it1 (outputImage, outputImage->GetLargestPossibleRegion()); boost::progress_display disp(outputImage->GetLargestPossibleRegion().GetNumberOfPixels()); while(!it1.IsAtEnd()) { if (this->GetAbortGenerateData()) { m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n"; return; } ++disp; unsigned long newTick = 50*disp.count()/disp.expected_count(); for (unsigned int tick = 0; tick<(newTick-lastTick); tick++) m_StatusText += "*"; lastTick = newTick; typename InputImageType::PixelType signal = it1.Get(); m_Parameters.m_NoiseModel->AddNoise(signal); it1.Set(signal); ++it1; } m_StatusText += "\n\n"; } this->SetNthOutput(0, outputImage); m_StatusText += "Finished simulation\n"; m_StatusText += "Simulation time: "+GetTime(); } template< class TPixelType > std::string AddArtifactsToDwiImageFilter< TPixelType >::GetTime() { unsigned long total = (double)(clock() - m_StartTime)/CLOCKS_PER_SEC; unsigned long hours = total/3600; unsigned long minutes = (total%3600)/60; unsigned long seconds = total%60; std::string out = ""; out.append(boost::lexical_cast(hours)); out.append(":"); out.append(boost::lexical_cast(minutes)); out.append(":"); out.append(boost::lexical_cast(seconds)); return out; } } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.cpp index 6888736d7e..77448fc7c4 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.cpp @@ -1,220 +1,220 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkKspaceImageFilter_txx #define __itkKspaceImageFilter_txx #include #include #include #include "itkKspaceImageFilter.h" #include #include #include #include #define _USE_MATH_DEFINES #include namespace itk { template< class TPixelType > KspaceImageFilter< TPixelType > ::KspaceImageFilter() : m_FrequencyMapSlice(NULL) , m_Z(0) , m_UseConstantRandSeed(false) , m_SpikesPerSlice(0) , m_IsBaseline(true) { m_DiffusionGradientDirection.Fill(0.0); m_RandGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); m_RandGen->SetSeed(); } template< class TPixelType > void KspaceImageFilter< TPixelType > ::BeforeThreadedGenerateData() { if (m_UseConstantRandSeed) // always generate the same random numbers? m_RandGen->SetSeed(0); else m_RandGen->SetSeed(); typename OutputImageType::Pointer outputImage = OutputImageType::New(); itk::ImageRegion<2> region; region.SetSize(0, m_OutSize[0]); region.SetSize(1, m_OutSize[1]); outputImage->SetLargestPossibleRegion( region ); outputImage->SetBufferedRegion( region ); outputImage->SetRequestedRegion( region ); outputImage->Allocate(); double gamma = 42576000; // Gyromagnetic ratio in Hz/T - if (m_Parameters.m_EddyStrength>0 && m_DiffusionGradientDirection.GetNorm()>0.001) + if ( m_Parameters.m_SignalGen.m_EddyStrength>0 && m_DiffusionGradientDirection.GetNorm()>0.001) { m_DiffusionGradientDirection.Normalize(); - m_DiffusionGradientDirection = m_DiffusionGradientDirection * m_Parameters.m_EddyStrength/1000 * gamma; + m_DiffusionGradientDirection = m_DiffusionGradientDirection * m_Parameters.m_SignalGen.m_EddyStrength/1000 * gamma; m_IsBaseline = false; } this->SetNthOutput(0, outputImage); m_Spike = vcl_complex(0,0); - m_Transform = m_Parameters.m_ImageDirection; + 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_ImageSpacing[j]; + m_Transform[i][j] *= m_Parameters.m_SignalGen.m_ImageSpacing[j]; } template< class TPixelType > void KspaceImageFilter< TPixelType > ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType) { typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); typedef ImageRegionConstIterator< InputImageType > InputIteratorType; double kxMax = outputImage->GetLargestPossibleRegion().GetSize(0); // k-space size in x-direction double kyMax = outputImage->GetLargestPossibleRegion().GetSize(1); // k-space size in y-direction double xMax = m_CompartmentImages.at(0)->GetLargestPossibleRegion().GetSize(0); // scanner coverage in x-direction double yMax = m_CompartmentImages.at(0)->GetLargestPossibleRegion().GetSize(1); // scanner coverage in y-direction double numPix = kxMax*kyMax; - double dt = m_Parameters.m_tLine/kxMax; - double fromMaxEcho = - m_Parameters.m_tLine*kyMax/2; + double dt = m_Parameters.m_SignalGen.m_tLine/kxMax; + double fromMaxEcho = - m_Parameters.m_SignalGen.m_tLine*kyMax/2; double upsampling = xMax/kxMax; // discrepany between k-space resolution and image resolution double yMaxFov = kyMax*upsampling; // actual FOV in y-direction (in x-direction xMax==FOV) int xRingingOffset = xMax-kxMax; int yRingingOffset = yMaxFov-kyMax; while( !oit.IsAtEnd() ) { itk::Index< 2 > kIdx; kIdx[0] = oit.GetIndex()[0]; kIdx[1] = oit.GetIndex()[1]; double t = fromMaxEcho + ((double)kIdx[1]*kxMax+(double)kIdx[0])*dt; // dephasing time // rearrange slice if( kIdx[0] < kxMax/2 ) kIdx[0] = kIdx[0] + kxMax/2; else kIdx[0] = kIdx[0] - kxMax/2; if( kIdx[1] < kyMax/2 ) kIdx[1] = kIdx[1] + kyMax/2; else kIdx[1] = kIdx[1] - kyMax/2; // calculate eddy current decay factors double eddyDecay = 0; - if (m_Parameters.m_EddyStrength>0) - eddyDecay = exp(-(m_Parameters.m_tEcho/2 + t)/m_Parameters.m_Tau) * t/1000; + if ( m_Parameters.m_SignalGen.m_EddyStrength>0) + eddyDecay = exp(-( m_Parameters.m_SignalGen.m_tEcho/2 + t)/ m_Parameters.m_SignalGen.m_Tau) * t/1000; // calcualte signal relaxation factors std::vector< double > relaxFactor; - if (m_Parameters.m_DoSimulateRelaxation) + if ( m_Parameters.m_SignalGen.m_DoSimulateRelaxation) for (unsigned int i=0; i=kxMax/2) kx += xRingingOffset; if (ky>=kyMax/2) ky += yRingingOffset; vcl_complex s(0,0); InputIteratorType it(m_CompartmentImages.at(0), m_CompartmentImages.at(0)->GetLargestPossibleRegion() ); while( !it.IsAtEnd() ) { double x = it.GetIndex()[0]-xMax/2; double y = it.GetIndex()[1]-yMax/2; vcl_complex f(0, 0); // sum compartment signals and simulate relaxation for (unsigned int i=0; i( m_CompartmentImages.at(i)->GetPixel(it.GetIndex()) * relaxFactor.at(i) * m_Parameters.m_SignalScale, 0); + if ( m_Parameters.m_SignalGen.m_DoSimulateRelaxation) + f += std::complex( 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_SignalScale ); + f += std::complex( m_CompartmentImages.at(i)->GetPixel(it.GetIndex()) * m_Parameters.m_SignalGen.m_SignalScale ); // simulate eddy currents and other distortions double omega_t = 0; - if ( m_Parameters.m_EddyStrength>0 && !m_IsBaseline) + if ( m_Parameters.m_SignalGen.m_EddyStrength>0 && !m_IsBaseline) { itk::Vector< double, 3 > pos; pos[0] = x; pos[1] = y; pos[2] = m_Z; pos = m_Transform*pos/1000; // vector from image center to current position (in meter) omega_t += (m_DiffusionGradientDirection[0]*pos[0]+m_DiffusionGradientDirection[1]*pos[1]+m_DiffusionGradientDirection[2]*pos[2])*eddyDecay; } if (m_FrequencyMapSlice.IsNotNull()) // simulate distortions omega_t += m_FrequencyMapSlice->GetPixel(it.GetIndex())*t/1000; if (y<-yMaxFov/2) y += yMaxFov; else if (y>=yMaxFov/2) y -= yMaxFov; // actual DFT term s += f * exp( std::complex(0, 2 * M_PI * (kx*x/xMax + ky*y/yMaxFov + omega_t )) ); ++it; } s /= numPix; if (m_SpikesPerSlice>0 && sqrt(s.imag()*s.imag()+s.real()*s.real()) > sqrt(m_Spike.imag()*m_Spike.imag()+m_Spike.real()*m_Spike.real()) ) m_Spike = s; outputImage->SetPixel(kIdx, s); ++oit; } } template< class TPixelType > void KspaceImageFilter< TPixelType > ::AfterThreadedGenerateData() { typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); double kxMax = outputImage->GetLargestPossibleRegion().GetSize(0); // k-space size in x-direction double kyMax = outputImage->GetLargestPossibleRegion().GetSize(1); // k-space size in y-direction - m_Spike *= m_Parameters.m_SpikeAmplitude; + m_Spike *= m_Parameters.m_SignalGen.m_SpikeAmplitude; itk::Index< 2 > spikeIdx; for (unsigned int i=0; iGetIntegerVariate()%(int)kxMax; spikeIdx[1] = m_RandGen->GetIntegerVariate()%(int)kyMax; outputImage->SetPixel(spikeIdx, m_Spike); } } } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp index f28b1cfe97..cf1803df39 100755 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp @@ -1,1300 +1,1300 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "itkTractsToDWIImageFilter.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include namespace itk { template< class PixelType > TractsToDWIImageFilter< PixelType >::TractsToDWIImageFilter() : m_FiberBundle(NULL) , m_StatusText("") , m_UseConstantRandSeed(false) , m_RandGen(itk::Statistics::MersenneTwisterRandomVariateGenerator::New()) { m_RandGen->SetSeed(); } template< class PixelType > TractsToDWIImageFilter< PixelType >::~TractsToDWIImageFilter() { } template< class PixelType > TractsToDWIImageFilter< PixelType >::DoubleDwiType::Pointer TractsToDWIImageFilter< PixelType >::DoKspaceStuff( std::vector< DoubleDwiType::Pointer >& images ) { int numFiberCompartments = m_Parameters.m_FiberModelList.size(); // create slice object ImageRegion<2> sliceRegion; sliceRegion.SetSize(0, m_UpsampledImageRegion.GetSize()[0]); sliceRegion.SetSize(1, m_UpsampledImageRegion.GetSize()[1]); Vector< double, 2 > sliceSpacing; sliceSpacing[0] = m_UpsampledSpacing[0]; sliceSpacing[1] = m_UpsampledSpacing[1]; // frequency map slice SliceType::Pointer fMapSlice = NULL; - if (m_Parameters.m_FrequencyMap.IsNotNull()) + if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull()) { fMapSlice = SliceType::New(); ImageRegion<2> region; region.SetSize(0, m_UpsampledImageRegion.GetSize()[0]); region.SetSize(1, m_UpsampledImageRegion.GetSize()[1]); fMapSlice->SetLargestPossibleRegion( region ); fMapSlice->SetBufferedRegion( region ); fMapSlice->SetRequestedRegion( region ); fMapSlice->Allocate(); fMapSlice->FillBuffer(0.0); } DoubleDwiType::Pointer newImage = DoubleDwiType::New(); - newImage->SetSpacing( m_Parameters.m_ImageSpacing ); - newImage->SetOrigin( m_Parameters.m_ImageOrigin ); - newImage->SetDirection( m_Parameters.m_ImageDirection ); - newImage->SetLargestPossibleRegion( m_Parameters.m_ImageRegion ); - newImage->SetBufferedRegion( m_Parameters.m_ImageRegion ); - newImage->SetRequestedRegion( m_Parameters.m_ImageRegion ); + newImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); + newImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin ); + newImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); + newImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_ImageRegion ); + newImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_ImageRegion ); + newImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_ImageRegion ); newImage->SetVectorLength( images.at(0)->GetVectorLength() ); newImage->Allocate(); std::vector< unsigned int > spikeVolume; - for (unsigned int i=0; iGetIntegerVariate()%images.at(0)->GetVectorLength()); std::sort (spikeVolume.begin(), spikeVolume.end()); std::reverse (spikeVolume.begin(), spikeVolume.end()); m_StatusText += "0% 10 20 30 40 50 60 70 80 90 100%\n"; m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*"; unsigned long lastTick = 0; boost::progress_display disp(2*images.at(0)->GetVectorLength()*images.at(0)->GetLargestPossibleRegion().GetSize(2)); for (unsigned int g=0; gGetVectorLength(); g++) { std::vector< unsigned int > spikeSlice; while (!spikeVolume.empty() && spikeVolume.back()==g) { spikeSlice.push_back(m_RandGen->GetIntegerVariate()%images.at(0)->GetLargestPossibleRegion().GetSize(2)); spikeVolume.pop_back(); } std::sort (spikeSlice.begin(), spikeSlice.end()); std::reverse (spikeSlice.begin(), spikeSlice.end()); for (unsigned int z=0; zGetLargestPossibleRegion().GetSize(2); z++) { std::vector< SliceType::Pointer > compartmentSlices; std::vector< double > t2Vector; for (unsigned int i=0; i* signalModel; if (iSetLargestPossibleRegion( sliceRegion ); slice->SetBufferedRegion( sliceRegion ); slice->SetRequestedRegion( sliceRegion ); slice->SetSpacing(sliceSpacing); slice->Allocate(); slice->FillBuffer(0.0); // extract slice from channel g for (unsigned int y=0; yGetLargestPossibleRegion().GetSize(1); y++) for (unsigned int x=0; xGetLargestPossibleRegion().GetSize(0); x++) { SliceType::IndexType index2D; index2D[0]=x; index2D[1]=y; DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z; slice->SetPixel(index2D, images.at(i)->GetPixel(index3D)[g]); if (fMapSlice.IsNotNull() && i==0) - fMapSlice->SetPixel(index2D, m_Parameters.m_FrequencyMap->GetPixel(index3D)); + fMapSlice->SetPixel(index2D, m_Parameters.m_SignalGen.m_FrequencyMap->GetPixel(index3D)); } compartmentSlices.push_back(slice); t2Vector.push_back(signalModel->GetT2()); } if (this->GetAbortGenerateData()) return NULL; // create k-sapce (inverse fourier transform slices) - itk::Size<2> outSize; outSize.SetElement(0, m_Parameters.m_ImageRegion.GetSize(0)); outSize.SetElement(1, m_Parameters.m_ImageRegion.GetSize(1)); + itk::Size<2> outSize; outSize.SetElement(0, m_Parameters.m_SignalGen.m_ImageRegion.GetSize(0)); outSize.SetElement(1, m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1)); itk::KspaceImageFilter< SliceType::PixelType >::Pointer idft = itk::KspaceImageFilter< SliceType::PixelType >::New(); idft->SetCompartmentImages(compartmentSlices); idft->SetT2(t2Vector); idft->SetUseConstantRandSeed(m_UseConstantRandSeed); idft->SetParameters(m_Parameters); idft->SetZ((double)z-(double)images.at(0)->GetLargestPossibleRegion().GetSize(2)/2.0); idft->SetDiffusionGradientDirection(m_Parameters.GetGradientDirection(g)); idft->SetFrequencyMapSlice(fMapSlice); idft->SetOutSize(outSize); int numSpikes = 0; while (!spikeSlice.empty() && spikeSlice.back()==z) { numSpikes++; spikeSlice.pop_back(); } idft->SetSpikesPerSlice(numSpikes); idft->Update(); ComplexSliceType::Pointer fSlice; fSlice = idft->GetOutput(); ++disp; unsigned long newTick = 50*disp.count()/disp.expected_count(); for (unsigned long tick = 0; tick<(newTick-lastTick); tick++) m_StatusText += "*"; lastTick = newTick; // fourier transform slice SliceType::Pointer newSlice; itk::DftImageFilter< SliceType::PixelType >::Pointer dft = itk::DftImageFilter< SliceType::PixelType >::New(); dft->SetInput(fSlice); dft->Update(); newSlice = dft->GetOutput(); // put slice back into channel g for (unsigned int y=0; yGetLargestPossibleRegion().GetSize(1); y++) for (unsigned int x=0; xGetLargestPossibleRegion().GetSize(0); x++) { DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z; SliceType::IndexType index2D; index2D[0]=x; index2D[1]=y; DoubleDwiType::PixelType pix3D = newImage->GetPixel(index3D); pix3D[g] = newSlice->GetPixel(index2D); newImage->SetPixel(index3D, pix3D); } ++disp; newTick = 50*disp.count()/disp.expected_count(); for (unsigned long tick = 0; tick<(newTick-lastTick); tick++) m_StatusText += "*"; lastTick = newTick; } } m_StatusText += "\n\n"; return newImage; } template< class PixelType > void TractsToDWIImageFilter< PixelType >::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!"); if (m_Parameters.m_NonFiberModelList.empty()) itkExceptionMacro("No diffusion model for non-fiber compartments defined!"); int baselineIndex = m_Parameters.GetFirstBaselineIndex(); if (baselineIndex<0) itkExceptionMacro("No baseline index found!"); - if (!m_Parameters.m_SimulateKspaceAcquisition) - m_Parameters.m_DoAddGibbsRinging = false; + if (!m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition) + m_Parameters.m_SignalGen.m_DoAddGibbsRinging = false; if (m_UseConstantRandSeed) // always generate the same random numbers? m_RandGen->SetSeed(0); else m_RandGen->SetSeed(); // initialize output dwi image - ImageRegion<3> croppedRegion = m_Parameters.m_ImageRegion; croppedRegion.SetSize(1, croppedRegion.GetSize(1)*m_Parameters.m_CroppingFactor); - itk::Point shiftedOrigin = m_Parameters.m_ImageOrigin; shiftedOrigin[1] += (m_Parameters.m_ImageRegion.GetSize(1)-croppedRegion.GetSize(1))*m_Parameters.m_ImageSpacing[1]/2; + ImageRegion<3> croppedRegion = m_Parameters.m_SignalGen.m_ImageRegion; croppedRegion.SetSize(1, 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)-croppedRegion.GetSize(1))*m_Parameters.m_SignalGen.m_ImageSpacing[1]/2; typename OutputImageType::Pointer outImage = OutputImageType::New(); - outImage->SetSpacing( m_Parameters.m_ImageSpacing ); + outImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); outImage->SetOrigin( shiftedOrigin ); - outImage->SetDirection( m_Parameters.m_ImageDirection ); + outImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); outImage->SetLargestPossibleRegion( croppedRegion ); outImage->SetBufferedRegion( croppedRegion ); outImage->SetRequestedRegion( croppedRegion ); outImage->SetVectorLength( m_Parameters.GetNumVolumes() ); outImage->Allocate(); typename OutputImageType::PixelType temp; temp.SetSize(m_Parameters.GetNumVolumes()); temp.Fill(0.0); outImage->FillBuffer(temp); // ADJUST GEOMETRY FOR FURTHER PROCESSING // is input slize size a power of two? - unsigned int x=m_Parameters.m_ImageRegion.GetSize(0); unsigned int y=m_Parameters.m_ImageRegion.GetSize(1); + 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_ImageRegion.SetSize(0, x+pad[0]); - m_Parameters.m_ImageRegion.SetSize(1, y+pad[1]); - if (m_Parameters.m_FrequencyMap.IsNotNull() && (pad[0]>0 || pad[1]>0)) + 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_FrequencyMap); + zeroPadder->SetInput(m_Parameters.m_SignalGen.m_FrequencyMap); zeroPadder->SetConstant(0); zeroPadder->SetPadUpperBound(pad); zeroPadder->Update(); - m_Parameters.m_FrequencyMap = zeroPadder->GetOutput(); + m_Parameters.m_SignalGen.m_FrequencyMap = zeroPadder->GetOutput(); } - if (m_Parameters.m_MaskImage.IsNotNull() && (pad[0]>0 || pad[1]>0)) + 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_MaskImage); + zeroPadder->SetInput(m_Parameters.m_SignalGen.m_MaskImage); zeroPadder->SetConstant(0); zeroPadder->SetPadUpperBound(pad); zeroPadder->Update(); - m_Parameters.m_MaskImage = zeroPadder->GetOutput(); + m_Parameters.m_SignalGen.m_MaskImage = zeroPadder->GetOutput(); } // Apply in-plane upsampling for Gibbs ringing artifact double upsampling = 1; - if (m_Parameters.m_DoAddGibbsRinging) + if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging) upsampling = 2; - m_UpsampledSpacing = m_Parameters.m_ImageSpacing; + m_UpsampledSpacing = m_Parameters.m_SignalGen.m_ImageSpacing; m_UpsampledSpacing[0] /= upsampling; m_UpsampledSpacing[1] /= upsampling; - m_UpsampledImageRegion = m_Parameters.m_ImageRegion; - m_UpsampledImageRegion.SetSize(0, m_Parameters.m_ImageRegion.GetSize()[0]*upsampling); - m_UpsampledImageRegion.SetSize(1, m_Parameters.m_ImageRegion.GetSize()[1]*upsampling); - m_UpsampledOrigin = m_Parameters.m_ImageOrigin; - m_UpsampledOrigin[0] -= m_Parameters.m_ImageSpacing[0]/2; m_UpsampledOrigin[0] += m_UpsampledSpacing[0]/2; - m_UpsampledOrigin[1] -= m_Parameters.m_ImageSpacing[1]/2; m_UpsampledOrigin[1] += m_UpsampledSpacing[1]/2; - m_UpsampledOrigin[2] -= m_Parameters.m_ImageSpacing[2]/2; m_UpsampledOrigin[2] += m_UpsampledSpacing[2]/2; + m_UpsampledImageRegion = m_Parameters.m_SignalGen.m_ImageRegion; + m_UpsampledImageRegion.SetSize(0, m_Parameters.m_SignalGen.m_ImageRegion.GetSize()[0]*upsampling); + m_UpsampledImageRegion.SetSize(1, m_Parameters.m_SignalGen.m_ImageRegion.GetSize()[1]*upsampling); + m_UpsampledOrigin = m_Parameters.m_SignalGen.m_ImageOrigin; + m_UpsampledOrigin[0] -= m_Parameters.m_SignalGen.m_ImageSpacing[0]/2; m_UpsampledOrigin[0] += m_UpsampledSpacing[0]/2; + m_UpsampledOrigin[1] -= m_Parameters.m_SignalGen.m_ImageSpacing[1]/2; m_UpsampledOrigin[1] += m_UpsampledSpacing[1]/2; + m_UpsampledOrigin[2] -= m_Parameters.m_SignalGen.m_ImageSpacing[2]/2; m_UpsampledOrigin[2] += m_UpsampledSpacing[2]/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_UpsampledSpacing ); doubleDwi->SetOrigin( m_UpsampledOrigin ); - doubleDwi->SetDirection( m_Parameters.m_ImageDirection ); + doubleDwi->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); doubleDwi->SetLargestPossibleRegion( m_UpsampledImageRegion ); doubleDwi->SetBufferedRegion( m_UpsampledImageRegion ); doubleDwi->SetRequestedRegion( m_UpsampledImageRegion ); doubleDwi->SetVectorLength( m_Parameters.GetNumVolumes() ); doubleDwi->Allocate(); DoubleDwiType::PixelType pix; pix.SetSize(m_Parameters.GetNumVolumes()); pix.Fill(0.0); doubleDwi->FillBuffer(pix); m_CompartmentImages.push_back(doubleDwi); } // initialize output volume fraction images m_VolumeFractions.clear(); for (int i=0; iSetSpacing( m_UpsampledSpacing ); doubleImg->SetOrigin( m_UpsampledOrigin ); - doubleImg->SetDirection( m_Parameters.m_ImageDirection ); + doubleImg->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); doubleImg->SetLargestPossibleRegion( m_UpsampledImageRegion ); doubleImg->SetBufferedRegion( m_UpsampledImageRegion ); doubleImg->SetRequestedRegion( m_UpsampledImageRegion ); doubleImg->Allocate(); doubleImg->FillBuffer(0); m_VolumeFractions.push_back(doubleImg); } // get volume fraction images ItkDoubleImgType::Pointer sumImage = ItkDoubleImgType::New(); bool foundVolumeFractionImage = false; for (int i=0; iGetVolumeFractionImage().IsNotNull()) { foundVolumeFractionImage = true; itk::ConstantPadImageFilter::Pointer zeroPadder = itk::ConstantPadImageFilter::New(); zeroPadder->SetInput(m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()); zeroPadder->SetConstant(0); zeroPadder->SetPadUpperBound(pad); zeroPadder->Update(); m_Parameters.m_NonFiberModelList[i]->SetVolumeFractionImage(zeroPadder->GetOutput()); sumImage->SetSpacing( m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetSpacing() ); sumImage->SetOrigin( m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetOrigin() ); sumImage->SetDirection( m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetDirection() ); sumImage->SetLargestPossibleRegion( m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetLargestPossibleRegion() ); sumImage->SetBufferedRegion( m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetLargestPossibleRegion() ); sumImage->SetRequestedRegion( m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetLargestPossibleRegion() ); sumImage->Allocate(); sumImage->FillBuffer(0); break; } } if (!foundVolumeFractionImage) { sumImage->SetSpacing( m_UpsampledSpacing ); sumImage->SetOrigin( m_UpsampledOrigin ); - sumImage->SetDirection( m_Parameters.m_ImageDirection ); + sumImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); sumImage->SetLargestPossibleRegion( m_UpsampledImageRegion ); sumImage->SetBufferedRegion( m_UpsampledImageRegion ); sumImage->SetRequestedRegion( m_UpsampledImageRegion ); sumImage->Allocate(); sumImage->FillBuffer(0.0); } for (int i=0; iGetVolumeFractionImage().IsNull()) { ItkDoubleImgType::Pointer doubleImg = ItkDoubleImgType::New(); doubleImg->SetSpacing( sumImage->GetSpacing() ); doubleImg->SetOrigin( sumImage->GetOrigin() ); doubleImg->SetDirection( sumImage->GetDirection() ); doubleImg->SetLargestPossibleRegion( sumImage->GetLargestPossibleRegion() ); doubleImg->SetBufferedRegion( sumImage->GetLargestPossibleRegion() ); doubleImg->SetRequestedRegion( sumImage->GetLargestPossibleRegion() ); doubleImg->Allocate(); doubleImg->FillBuffer(1.0/numNonFiberCompartments); m_Parameters.m_NonFiberModelList[i]->SetVolumeFractionImage(doubleImg); } ImageRegionIterator it(m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage(), m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetLargestPossibleRegion()); while(!it.IsAtEnd()) { sumImage->SetPixel(it.GetIndex(), sumImage->GetPixel(it.GetIndex())+it.Get()); ++it; } } for (int i=0; i it(m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage(), m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetLargestPossibleRegion()); while(!it.IsAtEnd()) { if (sumImage->GetPixel(it.GetIndex())>0) it.Set(it.Get()/sumImage->GetPixel(it.GetIndex())); ++it; } } // resample mask image and frequency map to fit upsampled geometry - if (m_Parameters.m_DoAddGibbsRinging) + if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging) { - if (m_Parameters.m_MaskImage.IsNotNull()) + 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_MaskImage); + 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_MaskImage); + resampler->SetOutputParametersFromImage(m_Parameters.m_SignalGen.m_MaskImage); resampler->SetSize(m_UpsampledImageRegion.GetSize()); resampler->SetOutputSpacing(m_UpsampledSpacing); resampler->SetOutputOrigin(m_UpsampledOrigin); itk::NearestNeighborInterpolateImageFunction::Pointer nn_interpolator = itk::NearestNeighborInterpolateImageFunction::New(); resampler->SetInterpolator(nn_interpolator); resampler->Update(); - m_Parameters.m_MaskImage = resampler->GetOutput(); + m_Parameters.m_SignalGen.m_MaskImage = resampler->GetOutput(); itk::ImageFileWriter::Pointer w = itk::ImageFileWriter::New(); w->SetFileName("/local/mask_ups.nrrd"); - w->SetInput(m_Parameters.m_MaskImage); + w->SetInput(m_Parameters.m_SignalGen.m_MaskImage); w->Update(); } // resample frequency map - if (m_Parameters.m_FrequencyMap.IsNotNull()) + if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull()) { itk::ResampleImageFilter::Pointer resampler = itk::ResampleImageFilter::New(); - resampler->SetInput(m_Parameters.m_FrequencyMap); - resampler->SetOutputParametersFromImage(m_Parameters.m_FrequencyMap); + resampler->SetInput(m_Parameters.m_SignalGen.m_FrequencyMap); + resampler->SetOutputParametersFromImage(m_Parameters.m_SignalGen.m_FrequencyMap); resampler->SetSize(m_UpsampledImageRegion.GetSize()); resampler->SetOutputSpacing(m_UpsampledSpacing); resampler->SetOutputOrigin(m_UpsampledOrigin); itk::NearestNeighborInterpolateImageFunction::Pointer nn_interpolator = itk::NearestNeighborInterpolateImageFunction::New(); resampler->SetInterpolator(nn_interpolator); resampler->Update(); - m_Parameters.m_FrequencyMap = resampler->GetOutput(); + m_Parameters.m_SignalGen.m_FrequencyMap = resampler->GetOutput(); } } // no input tissue mask is set -> create default m_MaskImageSet = true; - if (m_Parameters.m_MaskImage.IsNull()) + if (m_Parameters.m_SignalGen.m_MaskImage.IsNull()) { m_StatusText += "No tissue mask set\n"; MITK_INFO << "No tissue mask set"; - m_Parameters.m_MaskImage = ItkUcharImgType::New(); - m_Parameters.m_MaskImage->SetSpacing( m_UpsampledSpacing ); - m_Parameters.m_MaskImage->SetOrigin( m_UpsampledOrigin ); - m_Parameters.m_MaskImage->SetDirection( m_Parameters.m_ImageDirection ); - m_Parameters.m_MaskImage->SetLargestPossibleRegion( m_UpsampledImageRegion ); - m_Parameters.m_MaskImage->SetBufferedRegion( m_UpsampledImageRegion ); - m_Parameters.m_MaskImage->SetRequestedRegion( m_UpsampledImageRegion ); - m_Parameters.m_MaskImage->Allocate(); - m_Parameters.m_MaskImage->FillBuffer(1); + m_Parameters.m_SignalGen.m_MaskImage = ItkUcharImgType::New(); + m_Parameters.m_SignalGen.m_MaskImage->SetSpacing( m_UpsampledSpacing ); + m_Parameters.m_SignalGen.m_MaskImage->SetOrigin( m_UpsampledOrigin ); + m_Parameters.m_SignalGen.m_MaskImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); + m_Parameters.m_SignalGen.m_MaskImage->SetLargestPossibleRegion( m_UpsampledImageRegion ); + m_Parameters.m_SignalGen.m_MaskImage->SetBufferedRegion( m_UpsampledImageRegion ); + m_Parameters.m_SignalGen.m_MaskImage->SetRequestedRegion( m_UpsampledImageRegion ); + m_Parameters.m_SignalGen.m_MaskImage->Allocate(); + m_Parameters.m_SignalGen.m_MaskImage->FillBuffer(1); m_MaskImageSet = false; } else { m_StatusText += "Using tissue mask\n"; MITK_INFO << "Using tissue mask"; } - m_Parameters.m_ImageRegion = croppedRegion; - x=m_Parameters.m_ImageRegion.GetSize(0); y=m_Parameters.m_ImageRegion.GetSize(1); + m_Parameters.m_SignalGen.m_ImageRegion = croppedRegion; + x=m_Parameters.m_SignalGen.m_ImageRegion.GetSize(0); y=m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1); if ( x%2 == 1 ) - m_Parameters.m_ImageRegion.SetSize(0, x+1); + m_Parameters.m_SignalGen.m_ImageRegion.SetSize(0, x+1); if ( y%2 == 1 ) - m_Parameters.m_ImageRegion.SetSize(1, y+1); + m_Parameters.m_SignalGen.m_ImageRegion.SetSize(1, y+1); // resample fiber bundle for sufficient voxel coverage m_StatusText += "\n"+this->GetTime()+" > Resampling fibers ...\n"; double segmentVolume = 0.0001; float minSpacing = 1; if(m_UpsampledSpacing[0]GetDeepCopy(); double volumeAccuracy = 10; m_FiberBundleWorkingCopy->ResampleFibers(minSpacing/volumeAccuracy); - double mmRadius = m_Parameters.m_AxonRadius/1000; + double mmRadius = m_Parameters.m_SignalGen.m_AxonRadius/1000; if (mmRadius>0) segmentVolume = M_PI*mmRadius*mmRadius*minSpacing/volumeAccuracy; double maxVolume = 0; m_VoxelVolume = m_UpsampledSpacing[0]*m_UpsampledSpacing[1]*m_UpsampledSpacing[2]; - if (m_Parameters.m_DoAddMotion) + if (m_Parameters.m_SignalGen.m_DoAddMotion) { std::string fileName = "fiberfox_motion_0.log"; std::string filePath = mitk::IOUtil::GetTempPath(); - if (m_Parameters.m_OutputPath.size()>0) - filePath = m_Parameters.m_OutputPath; + 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_DoRandomizeMotion) + if (m_Parameters.m_SignalGen.m_DoRandomizeMotion) { m_StatusText += "Adding random motion artifacts:\n"; - m_StatusText += "Maximum rotation: +/-" + boost::lexical_cast(m_Parameters.m_Rotation) + "°\n"; - m_StatusText += "Maximum translation: +/-" + boost::lexical_cast(m_Parameters.m_Translation) + "mm\n"; + 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_Rotation) + "°\n"; - m_StatusText += "Maximum translation: " + boost::lexical_cast(m_Parameters.m_Translation) + "mm\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_Rotation; - MITK_INFO << "Maxmimum translation: " << m_Parameters.m_Translation; + MITK_INFO << "Maximum rotation: " << m_Parameters.m_SignalGen.m_Rotation; + MITK_INFO << "Maxmimum translation: " << m_Parameters.m_SignalGen.m_Translation; } maxVolume = 0; m_StatusText += "\n"+this->GetTime()+" > Generating " + boost::lexical_cast(numFiberCompartments+numNonFiberCompartments) + "-compartment diffusion-weighted signal.\n"; int numFibers = m_FiberBundle->GetNumFibers(); boost::progress_display disp(numFibers*m_Parameters.GetNumVolumes()); // get transform for motion artifacts m_FiberBundleTransformed = m_FiberBundleWorkingCopy; - m_Rotation = m_Parameters.m_Rotation/m_Parameters.GetNumVolumes(); - m_Translation = m_Parameters.m_Translation/m_Parameters.GetNumVolumes(); + m_Rotation = m_Parameters.m_SignalGen.m_Rotation/m_Parameters.GetNumVolumes(); + m_Translation = m_Parameters.m_SignalGen.m_Translation/m_Parameters.GetNumVolumes(); // creat image to hold transformed mask (motion artifact) m_MaskImage = ItkUcharImgType::New(); itk::ImageDuplicator::Pointer duplicator = itk::ImageDuplicator::New(); - duplicator->SetInputImage(m_Parameters.m_MaskImage); + duplicator->SetInputImage(m_Parameters.m_SignalGen.m_MaskImage); duplicator->Update(); m_MaskImage = duplicator->GetOutput(); // second upsampling needed for motion artifacts ImageRegion<3> upsampledImageRegion = m_UpsampledImageRegion; DoubleVectorType upsampledSpacing = m_UpsampledSpacing; upsampledSpacing[0] /= 4; upsampledSpacing[1] /= 4; upsampledSpacing[2] /= 4; upsampledImageRegion.SetSize(0, m_UpsampledImageRegion.GetSize()[0]*4); upsampledImageRegion.SetSize(1, m_UpsampledImageRegion.GetSize()[1]*4); upsampledImageRegion.SetSize(2, m_UpsampledImageRegion.GetSize()[2]*4); itk::Point upsampledOrigin = m_UpsampledOrigin; upsampledOrigin[0] -= m_UpsampledSpacing[0]/2; upsampledOrigin[0] += upsampledSpacing[0]/2; upsampledOrigin[1] -= m_UpsampledSpacing[1]/2; upsampledOrigin[1] += upsampledSpacing[1]/2; upsampledOrigin[2] -= m_UpsampledSpacing[2]/2; upsampledOrigin[2] += upsampledSpacing[2]/2; m_UpsampledMaskImage = ItkUcharImgType::New(); itk::ResampleImageFilter::Pointer upsampler = itk::ResampleImageFilter::New(); - upsampler->SetInput(m_Parameters.m_MaskImage); - upsampler->SetOutputParametersFromImage(m_Parameters.m_MaskImage); + 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(); unsigned long lastTick = 0; int signalModelSeed = m_RandGen->GetIntegerVariate(); - switch (m_Parameters.m_DiffusionDirectionMode) + switch (m_Parameters.m_SignalGen.m_DiffusionDirectionMode) { - case(FiberfoxParameters<>::FIBER_TANGENT_DIRECTIONS): // use fiber tangent directions to determine diffusion direction + 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); ItkDoubleImgType::Pointer intraAxonalVolumeImage = ItkDoubleImgType::New(); intraAxonalVolumeImage->SetSpacing( m_UpsampledSpacing ); intraAxonalVolumeImage->SetOrigin( m_UpsampledOrigin ); - intraAxonalVolumeImage->SetDirection( m_Parameters.m_ImageDirection ); + intraAxonalVolumeImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); intraAxonalVolumeImage->SetLargestPossibleRegion( m_UpsampledImageRegion ); intraAxonalVolumeImage->SetBufferedRegion( m_UpsampledImageRegion ); intraAxonalVolumeImage->SetRequestedRegion( m_UpsampledImageRegion ); intraAxonalVolumeImage->Allocate(); intraAxonalVolumeImage->FillBuffer(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; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (numPoints<2) continue; for( int j=0; jGetAbortGenerateData()) { m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n"; return; } double* temp = points->GetPoint(j); itk::Point vertex = GetItkPoint(temp); itk::Vector v = GetItkVector(temp); itk::Vector dir(3); if (jGetPoint(j+1))-v; else dir = v-GetItkVector(points->GetPoint(j-1)); if (dir.GetSquaredNorm()<0.0001 || dir[0]!=dir[0] || dir[1]!=dir[1] || dir[2]!=dir[2]) continue; itk::Index<3> idx; itk::ContinuousIndex contIndex; m_MaskImage->TransformPhysicalPointToIndex(vertex, idx); m_MaskImage->TransformPhysicalPointToContinuousIndex(vertex, contIndex); if (!m_MaskImage->GetLargestPossibleRegion().IsInside(idx) || m_MaskImage->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] += segmentVolume*m_Parameters.m_FiberModelList[k]->SimulateMeasurement(g); m_CompartmentImages.at(k)->SetPixel(idx, pix); } // update fiber volume image double vol = intraAxonalVolumeImage->GetPixel(idx) + segmentVolume; intraAxonalVolumeImage->SetPixel(idx, vol); if (g==0 && vol>maxVolume) 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_MaskImage, m_MaskImage->GetLargestPossibleRegion()); double fact = 1; - if (m_Parameters.m_AxonRadius<0.0001 || maxVolume>m_VoxelVolume) + if (m_Parameters.m_SignalGen.m_AxonRadius<0.0001 || maxVolume>m_VoxelVolume) fact = m_VoxelVolume/maxVolume; while(!it3.IsAtEnd()) { if (it3.Get()>0) { DoubleDwiType::IndexType index = it3.GetIndex(); // adjust intra-axonal signal to abtain an only-fiber voxel if (fabs(fact-1.0)>0.0001) for (int i=0; iGetPixel(index); pix[g] *= fact; m_CompartmentImages.at(i)->SetPixel(index, pix); } // simulate other compartments SimulateNonFiberSignal(index, intraAxonalVolumeImage->GetPixel(index)*fact, g); } ++it3; } // move fibers SimulateMotion(g); } break; } - case (FiberfoxParameters<>::MAIN_FIBER_DIRECTIONS): // use main fiber directions to determine voxel-wise diffusion directions + 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_MaskImage); - fOdfFilter->SetAngularThreshold(cos(m_Parameters.m_FiberSeparationThreshold*M_PI/180.0)); + 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_UpsampledSpacing ); intraAxonalVolumeImage->SetOrigin( m_UpsampledOrigin ); - intraAxonalVolumeImage->SetDirection( m_Parameters.m_ImageDirection ); + intraAxonalVolumeImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); intraAxonalVolumeImage->SetLargestPossibleRegion( m_UpsampledImageRegion ); intraAxonalVolumeImage->SetBufferedRegion( m_UpsampledImageRegion ); intraAxonalVolumeImage->SetRequestedRegion( m_UpsampledImageRegion ); 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_MaskImage->GetLargestPossibleRegion().GetNumberOfPixels()*m_Parameters.GetNumVolumes()); for (unsigned int g=0; gSetSeed(signalModelSeed); for (int i=0; iSetSeed(signalModelSeed); - if (m_Parameters.m_DoAddMotion && g>0) // if fibers have moved we need a new TDI and new directions + 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_MaskImage); fOdfFilter->Update(); directionImageContainer = fOdfFilter->GetDirectionImageContainer(); tdiFilter->SetFiberBundle(m_FiberBundleTransformed); tdiFilter->Update(); intraAxonalVolumeImage = tdiFilter->GetOutput(); } ImageRegionIterator< ItkUcharImgType > it(m_MaskImage, m_MaskImage->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 SimulateNonFiberSignal(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 (FiberfoxParameters<>::RANDOM_DIRECTIONS): + case (SignalGenerationParameters::RANDOM_DIRECTIONS): { ItkUcharImgType::Pointer numDirectionsImage = ItkUcharImgType::New(); numDirectionsImage->SetSpacing( m_UpsampledSpacing ); numDirectionsImage->SetOrigin( m_UpsampledOrigin ); - numDirectionsImage->SetDirection( m_Parameters.m_ImageDirection ); + numDirectionsImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); numDirectionsImage->SetLargestPossibleRegion( m_UpsampledImageRegion ); numDirectionsImage->SetBufferedRegion( m_UpsampledImageRegion ); numDirectionsImage->SetRequestedRegion( m_UpsampledImageRegion ); numDirectionsImage->Allocate(); numDirectionsImage->FillBuffer(0); - double sepAngle = cos(m_Parameters.m_FiberSeparationThreshold*M_PI/180.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_MaskImage->GetLargestPossibleRegion().GetNumberOfPixels()); ImageRegionIterator it(m_MaskImage, m_MaskImage->GetLargestPossibleRegion()); while(!it.IsAtEnd()) { ++disp; unsigned long newTick = 50*disp.count()/disp.expected_count(); for (unsigned int tick = 0; tick<(newTick-lastTick); tick++) m_StatusText += "*"; lastTick = newTick; if (this->GetAbortGenerateData()) { m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n"; return; } if (it.Get()>0) { int numFibs = m_RandGen->GetIntegerVariate(2)+1; DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(it.GetIndex()); double volume = m_RandGen->GetVariateWithClosedRange(0.3); // double sum = 0; std::vector< double > fractions; for (int i=0; iGetVariateWithClosedRange(0.5)); // sum += fractions.at(i); } // for (int i=0; i > directions; for (int i=0; iGetVariateWithClosedRange(2)-1.0; fib[1] = m_RandGen->GetVariateWithClosedRange(2)-1.0; fib[2] = m_RandGen->GetVariateWithClosedRange(2)-1.0; fib.Normalize(); double min = 0; for (unsigned int d=0; dmin) min = angle; } if (minSetFiberDirection(fib); pix += m_Parameters.m_FiberModelList.at(0)->SimulateMeasurement()*fractions[i]; directions.push_back(fib); } else i--; } pix *= (1-volume); m_CompartmentImages.at(0)->SetPixel(it.GetIndex(), pix); // CSF/GM { pix += volume*m_Parameters.m_NonFiberModelList.at(0)->SimulateMeasurement(); } numDirectionsImage->SetPixel(it.GetIndex(), numFibs); } ++it; } itk::ImageFileWriter< ItkUcharImgType >::Pointer wr = itk::ImageFileWriter< ItkUcharImgType >::New(); wr->SetInput(numDirectionsImage); wr->SetFileName(mitk::IOUtil::GetTempPath()+"/NumDirections_RandomDirections.nrrd"); wr->Update(); } } if (m_Logfile.is_open()) { m_Logfile << "DONE"; m_Logfile.close(); } m_StatusText += "\n\n"; if (this->GetAbortGenerateData()) { m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n"; return; } DoubleDwiType::Pointer doubleOutImage; - if ( m_Parameters.m_SimulateKspaceAcquisition ) // do k-space stuff + if ( m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition ) // do k-space stuff { m_StatusText += this->GetTime()+" > Adjusting complex signal\n"; MITK_INFO << "Adjusting complex signal:"; - if (m_Parameters.m_DoSimulateRelaxation) + if (m_Parameters.m_SignalGen.m_DoSimulateRelaxation) m_StatusText += "Simulating signal relaxation\n"; - if (m_Parameters.m_FrequencyMap.IsNotNull()) + if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull()) m_StatusText += "Simulating distortions\n"; - if (m_Parameters.m_DoAddGibbsRinging) + if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging) m_StatusText += "Simulating ringing artifacts\n"; - if (m_Parameters.m_EddyStrength>0) + if (m_Parameters.m_SignalGen.m_EddyStrength>0) m_StatusText += "Simulating eddy currents\n"; - if (m_Parameters.m_Spikes>0) + if (m_Parameters.m_SignalGen.m_Spikes>0) m_StatusText += "Simulating spikes\n"; - if (m_Parameters.m_CroppingFactor<1.0) + if (m_Parameters.m_SignalGen.m_CroppingFactor<1.0) m_StatusText += "Simulating aliasing artifacts\n"; - if (m_Parameters.m_KspaceLineOffset>0) + if (m_Parameters.m_SignalGen.m_KspaceLineOffset>0) m_StatusText += "Simulating ghosts\n"; doubleOutImage = DoKspaceStuff(m_CompartmentImages); - m_Parameters.m_SignalScale = 1; // already scaled in DoKspaceStuff + m_Parameters.m_SignalGen.m_SignalScale = 1; // already scaled in DoKspaceStuff } else // don't do k-space stuff, just sum compartments { m_StatusText += this->GetTime()+" > Summing compartments\n"; MITK_INFO << "Summing compartments"; doubleOutImage = m_CompartmentImages.at(0); for (unsigned int i=1; i::Pointer adder = itk::AddImageFilter< DoubleDwiType, DoubleDwiType, DoubleDwiType>::New(); adder->SetInput1(doubleOutImage); adder->SetInput2(m_CompartmentImages.at(i)); adder->Update(); doubleOutImage = adder->GetOutput(); } } if (this->GetAbortGenerateData()) { m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n"; return; } m_StatusText += this->GetTime()+" > Finalizing image\n"; MITK_INFO << "Finalizing image"; - if (m_Parameters.m_SignalScale>1) + if (m_Parameters.m_SignalGen.m_SignalScale>1) m_StatusText += " Scaling signal\n"; if (m_Parameters.m_NoiseModel!=NULL) m_StatusText += " Adding noise\n"; unsigned int window = 0; unsigned int min = itk::NumericTraits::max(); ImageRegionIterator it4 (outImage, outImage->GetLargestPossibleRegion()); DoubleDwiType::PixelType signal; signal.SetSize(m_Parameters.GetNumVolumes()); boost::progress_display disp2(outImage->GetLargestPossibleRegion().GetNumberOfPixels()); m_StatusText += "0% 10 20 30 40 50 60 70 80 90 100%\n"; m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*"; lastTick = 0; while(!it4.IsAtEnd()) { if (this->GetAbortGenerateData()) { m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n"; return; } ++disp2; unsigned long newTick = 50*disp2.count()/disp2.expected_count(); for (unsigned long tick = 0; tick<(newTick-lastTick); tick++) m_StatusText += "*"; lastTick = newTick; typename OutputImageType::IndexType index = it4.GetIndex(); - signal = doubleOutImage->GetPixel(index)*m_Parameters.m_SignalScale; + signal = doubleOutImage->GetPixel(index)*m_Parameters.m_SignalGen.m_SignalScale; if (m_Parameters.m_NoiseModel!=NULL) m_Parameters.m_NoiseModel->AddNoise(signal); for (unsigned int i=0; i0) signal[i] = floor(signal[i]+0.5); else signal[i] = ceil(signal[i]-0.5); if ( (!m_Parameters.IsBaselineIndex(i) || signal.Size()==1) && signal[i]>window) window = signal[i]; if ( (!m_Parameters.IsBaselineIndex(i) || signal.Size()==1) && signal[i]SetNthOutput(0, outImage); 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_DoAddMotion && gGetDeepCopy(); - m_Rotation[0] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_Rotation[0]*2)-m_Parameters.m_Rotation[0]; - m_Rotation[1] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_Rotation[1]*2)-m_Parameters.m_Rotation[1]; - m_Rotation[2] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_Rotation[2]*2)-m_Parameters.m_Rotation[2]; - m_Translation[0] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_Translation[0]*2)-m_Parameters.m_Translation[0]; - m_Translation[1] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_Translation[1]*2)-m_Parameters.m_Translation[1]; - m_Translation[2] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_Translation[2]*2)-m_Parameters.m_Translation[2]; + 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_MaskImage->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_DoRandomizeMotion) + 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_MaskImage->TransformPhysicalPointToIndex(point, index); if (m_MaskImage->GetLargestPossibleRegion().IsInside(index)) m_MaskImage->SetPixel(index,100); ++maskIt; } } // rotate fibers if (m_Logfile.is_open()) { m_Logfile << g+1 << " rotation: " << m_Rotation[0] << "," << m_Rotation[1] << "," << m_Rotation[2] << ";"; m_Logfile << " translation: " << m_Translation[0] << "," << m_Translation[1] << "," << m_Translation[2] << "\n"; } m_FiberBundleTransformed->TransformFibers(m_Rotation[0],m_Rotation[1],m_Rotation[2],m_Translation[0],m_Translation[1],m_Translation[2]); } } template< class PixelType > void TractsToDWIImageFilter< PixelType >::SimulateNonFiberSignal(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_DoDisablePartialVolume) // only fiber in voxel + 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); 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 { m_VolumeFractions.at(0)->SetPixel(index, intraAxonalVolume/m_VoxelVolume); itk::Point point; m_MaskImage->TransformIndexToPhysicalPoint(index, point); - if (m_Parameters.m_DoAddMotion) + if (m_Parameters.m_SignalGen.m_DoAddMotion) { - if (m_Parameters.m_DoRandomizeMotion && g>0) + 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_DoDisablePartialVolume) + if (m_Parameters.m_SignalGen.m_DoDisablePartialVolume) { int maxVolumeIndex = 0; double maxWeight = 0; for (int i=0; i1) { 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); m_VolumeFractions.at(maxVolumeIndex+numFiberCompartments)->SetPixel(index, 1); } else { double extraAxonalVolume = m_VoxelVolume-intraAxonalVolume; // non-fiber volume double interAxonalVolume = 0; if (numFiberCompartments>1) interAxonalVolume = extraAxonalVolume * intraAxonalVolume/m_VoxelVolume; // inter-axonal fraction of non fiber compartment scales linearly with f double other = extraAxonalVolume - interAxonalVolume; // rest of compartment double singleinter = interAxonalVolume/(numFiberCompartments-1); // 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 (g>=0) pix[g] *= singleinter; else pix *= singleinter; m_CompartmentImages.at(i)->SetPixel(index, pix); m_VolumeFractions.at(i)->SetPixel(index, singleinter/m_VoxelVolume); } for (int i=0; i1) { 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); } DoubleDwiType::Pointer doubleDwi = m_CompartmentImages.at(i+numFiberCompartments); DoubleDwiType::PixelType pix = doubleDwi->GetPixel(index); if (g>=0) pix[g] += m_Parameters.m_NonFiberModelList[i]->SimulateMeasurement(g)*other*weight; else pix += m_Parameters.m_NonFiberModelList[i]->SimulateMeasurement()*other*weight; doubleDwi->SetPixel(index, pix); m_VolumeFractions.at(i+numFiberCompartments)->SetPixel(index, other/m_VoxelVolume*weight); } } } } template< class PixelType > itk::Point TractsToDWIImageFilter< PixelType >::GetItkPoint(double point[3]) { itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; return itkPoint; } template< class PixelType > itk::Vector TractsToDWIImageFilter< PixelType >::GetItkVector(double point[3]) { itk::Vector itkVector; itkVector[0] = point[0]; itkVector[1] = point[1]; itkVector[2] = point[2]; return itkVector; } template< class PixelType > vnl_vector_fixed TractsToDWIImageFilter< PixelType >::GetVnlVector(double point[3]) { vnl_vector_fixed vnlVector; vnlVector[0] = point[0]; vnlVector[1] = point[1]; vnlVector[2] = point[2]; return vnlVector; } template< class PixelType > vnl_vector_fixed TractsToDWIImageFilter< PixelType >::GetVnlVector(Vector& vector) { vnl_vector_fixed vnlVector; vnlVector[0] = vector[0]; vnlVector[1] = vector[1]; vnlVector[2] = vector[2]; return vnlVector; } template< class PixelType > double TractsToDWIImageFilter< PixelType >::RoundToNearest(double num) { return (num > 0.0) ? floor(num + 0.5) : ceil(num - 0.5); } template< class PixelType > std::string TractsToDWIImageFilter< PixelType >::GetTime() { m_TimeProbe.Stop(); unsigned long total = RoundToNearest(m_TimeProbe.GetTotal()); unsigned long hours = total/3600; unsigned long minutes = (total%3600)/60; unsigned long seconds = total%60; std::string out = ""; out.append(boost::lexical_cast(hours)); out.append(":"); out.append(boost::lexical_cast(minutes)); out.append(":"); out.append(boost::lexical_cast(seconds)); m_TimeProbe.Start(); return out; } } diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.cpp index 76873960c7..a12c14b88a 100644 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.cpp +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.cpp @@ -1,648 +1,648 @@ /*=================================================================== 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 template< class ScalarType > mitk::FiberfoxParameters< ScalarType >::FiberfoxParameters() - : m_SignalScale(100) - , m_tEcho(100) - , m_tLine(1) - , m_tInhom(50) - , m_Bvalue(1000) - , m_SimulateKspaceAcquisition(false) - , m_AxonRadius(0) - , m_DiffusionDirectionMode(FIBER_TANGENT_DIRECTIONS) - , m_FiberSeparationThreshold(30) - , m_Spikes(0) - , m_SpikeAmplitude(1) - , m_KspaceLineOffset(0) - , m_EddyStrength(0) - , m_Tau(70) - , m_CroppingFactor(1) - , m_DoAddNoise(false) - , m_DoAddGhosts(false) - , m_DoAddAliasing(false) - , m_DoAddSpikes(false) - , m_DoAddEddyCurrents(false) - , m_DoAddGibbsRinging(false) - , m_DoSimulateRelaxation(true) - , m_DoAddDistortions(false) - , m_DoDisablePartialVolume(false) - , m_DoAddMotion(false) - , m_DoRandomizeMotion(true) - , m_NoiseModel(NULL) - , m_FrequencyMap(NULL) - , m_MaskImage(NULL) - , m_ResultNode(mitk::DataNode::New()) - , m_ParentNode(NULL) - , m_SignalModelString("") - , m_ArtifactModelString("") - , m_OutputPath("") + : m_NoiseModel(NULL) , m_NumGradients(6) , m_NumBaseline(1) { - m_FiberGenerationParameters.m_RealTimeFibers = true; - m_FiberGenerationParameters.m_AdvancedOptions = false; - m_FiberGenerationParameters.m_Distribution = 0; - m_FiberGenerationParameters.m_Variance = 100; - m_FiberGenerationParameters.m_FiberDensity = 100; - m_FiberGenerationParameters.m_Sampling = 1; - m_FiberGenerationParameters.m_Tension = 0; - m_FiberGenerationParameters.m_Continuity = 0; - m_FiberGenerationParameters.m_Bias = 0; - m_FiberGenerationParameters.m_ConstantRadius = false; - m_FiberGenerationParameters.m_Rotation[0] = 0; - m_FiberGenerationParameters.m_Rotation[1] = 0; - m_FiberGenerationParameters.m_Rotation[2] = 0; - m_FiberGenerationParameters.m_Translation[0] = 0; - m_FiberGenerationParameters.m_Translation[1] = 0; - m_FiberGenerationParameters.m_Translation[2] = 0; - m_FiberGenerationParameters.m_Scale[0] = 1; - m_FiberGenerationParameters.m_Scale[1] = 1; - m_FiberGenerationParameters.m_Scale[2] = 1; - m_FiberGenerationParameters.m_IncludeFiducials = true; - - m_ImageDirection.SetIdentity(); - m_ImageOrigin.Fill(0.0); - m_ImageRegion.SetSize(0, 11); - m_ImageRegion.SetSize(1, 11); - m_ImageRegion.SetSize(2, 3); - m_ImageSpacing.Fill(2.0); - - m_Translation.Fill(0.0); - m_Rotation.Fill(0.0); + m_FiberGen.m_RealTimeFibers = true; + m_FiberGen.m_AdvancedOptions = false; + m_FiberGen.m_Distribution = 0; + m_FiberGen.m_Variance = 100; + m_FiberGen.m_FiberDensity = 100; + m_FiberGen.m_Sampling = 1; + m_FiberGen.m_Tension = 0; + m_FiberGen.m_Continuity = 0; + m_FiberGen.m_Bias = 0; + m_FiberGen.m_ConstantRadius = false; + m_FiberGen.m_Rotation[0] = 0; + m_FiberGen.m_Rotation[1] = 0; + m_FiberGen.m_Rotation[2] = 0; + m_FiberGen.m_Translation[0] = 0; + m_FiberGen.m_Translation[1] = 0; + m_FiberGen.m_Translation[2] = 0; + m_FiberGen.m_Scale[0] = 1; + m_FiberGen.m_Scale[1] = 1; + m_FiberGen.m_Scale[2] = 1; + m_FiberGen.m_IncludeFiducials = true; + + m_SignalGen.m_ImageDirection.SetIdentity(); + m_SignalGen.m_ImageOrigin.Fill(0.0); + m_SignalGen.m_ImageRegion.SetSize(0, 11); + m_SignalGen.m_ImageRegion.SetSize(1, 11); + m_SignalGen.m_ImageRegion.SetSize(2, 3); + m_SignalGen.m_ImageSpacing.Fill(2.0); + m_SignalGen.m_Translation.Fill(0.0); + m_SignalGen.m_Rotation.Fill(0.0); + m_SignalGen.m_FrequencyMap = NULL; + m_SignalGen.m_MaskImage = NULL; + m_SignalGen.m_SignalScale = 100; + m_SignalGen.m_tEcho = 100; + m_SignalGen.m_tLine = 1; + m_SignalGen.m_tInhom = 50; + m_SignalGen.m_Bvalue = 1000; + m_SignalGen.m_SimulateKspaceAcquisition = false; + m_SignalGen.m_AxonRadius = 0; + m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::FIBER_TANGENT_DIRECTIONS; + m_SignalGen.m_FiberSeparationThreshold = 30; + m_SignalGen.m_Spikes = 0; + m_SignalGen.m_SpikeAmplitude = 1; + m_SignalGen.m_KspaceLineOffset = 0; + m_SignalGen.m_EddyStrength = 0; + m_SignalGen.m_Tau = 70; + m_SignalGen.m_CroppingFactor = 1; + m_SignalGen.m_DoAddNoise = false; + m_SignalGen.m_DoAddGhosts = false; + m_SignalGen.m_DoAddAliasing = false; + m_SignalGen.m_DoAddSpikes = false; + m_SignalGen.m_DoAddEddyCurrents = false; + m_SignalGen.m_DoAddGibbsRinging = false; + m_SignalGen.m_DoSimulateRelaxation = true; + m_SignalGen.m_DoAddDistortions = false; + m_SignalGen.m_DoDisablePartialVolume = false; + m_SignalGen.m_DoAddMotion = false; + m_SignalGen.m_DoRandomizeMotion = true; + + m_Misc.m_ResultNode = mitk::DataNode::New(); + m_Misc.m_ParentNode = NULL; + m_Misc.m_SignalModelString = ""; + m_Misc.m_ArtifactModelString = ""; + m_Misc.m_OutputPath = ""; GenerateGradientHalfShell(); } template< class ScalarType > mitk::FiberfoxParameters< ScalarType >::~FiberfoxParameters() { // if (m_NoiseModel!=NULL) // delete m_NoiseModel; } template< class ScalarType > void mitk::FiberfoxParameters< ScalarType >::GenerateGradientHalfShell() { int NPoints = 2*m_NumGradients; m_GradientDirections.clear(); m_NumBaseline = NPoints/20; if (m_NumBaseline==0) m_NumBaseline=1; GradientType g; g.Fill(0.0); for (unsigned int i=0; i theta; theta.set_size(NPoints); vnl_vector phi; phi.set_size(NPoints); double C = sqrt(4*M_PI); phi(0) = 0.0; phi(NPoints-1) = 0.0; for(int i=0; i0 && i std::vector< int > mitk::FiberfoxParameters< ScalarType >::GetBaselineIndices() { std::vector< int > result; for( unsigned int i=0; im_GradientDirections.size(); i++) if (m_GradientDirections.at(i).GetNorm()<0.0001) result.push_back(i); return result; } template< class ScalarType > unsigned int mitk::FiberfoxParameters< ScalarType >::GetFirstBaselineIndex() { for( unsigned int i=0; im_GradientDirections.size(); i++) if (m_GradientDirections.at(i).GetNorm()<0.0001) return i; return -1; } template< class ScalarType > bool mitk::FiberfoxParameters< ScalarType >::IsBaselineIndex(unsigned int idx) { if (m_GradientDirections.size()>idx && m_GradientDirections.at(idx).GetNorm()<0.0001) return true; return false; } template< class ScalarType > unsigned int mitk::FiberfoxParameters< ScalarType >::GetNumWeightedVolumes() { return m_NumGradients; } template< class ScalarType > unsigned int mitk::FiberfoxParameters< ScalarType >::GetNumBaselineVolumes() { return m_NumBaseline; } template< class ScalarType > unsigned int mitk::FiberfoxParameters< ScalarType >::GetNumVolumes() { return m_GradientDirections.size(); } template< class ScalarType > typename mitk::FiberfoxParameters< ScalarType >::GradientListType mitk::FiberfoxParameters< ScalarType >::GetGradientDirections() { return m_GradientDirections; } template< class ScalarType > typename mitk::FiberfoxParameters< ScalarType >::GradientType mitk::FiberfoxParameters< ScalarType >::GetGradientDirection(unsigned int i) { return m_GradientDirections.at(i); } template< class ScalarType > void mitk::FiberfoxParameters< ScalarType >::SetNumWeightedGradients(int numGradients) { m_NumGradients = numGradients; GenerateGradientHalfShell(); } template< class ScalarType > void mitk::FiberfoxParameters< ScalarType >::SetGradienDirections(GradientListType gradientList) { m_GradientDirections = gradientList; m_NumGradients = 0; m_NumBaseline = 0; for( unsigned int i=0; im_GradientDirections.size(); i++) { if (m_GradientDirections.at(i).GetNorm()>0.0001) m_NumGradients++; else m_NumBaseline++; } } template< class ScalarType > void mitk::FiberfoxParameters< ScalarType >::SetGradienDirections(mitk::DiffusionImage::GradientDirectionContainerType::Pointer gradientList) { m_NumGradients = 0; m_NumBaseline = 0; m_GradientDirections.clear(); for( unsigned int i=0; iSize(); i++) { GradientType g; g[0] = gradientList->at(i)[0]; g[1] = gradientList->at(i)[1]; g[2] = gradientList->at(i)[2]; m_GradientDirections.push_back(g); if (m_GradientDirections.at(i).GetNorm()>0.0001) m_NumGradients++; else m_NumBaseline++; } } template< class ScalarType > void mitk::FiberfoxParameters< ScalarType >::SaveParameters(string filename) { if(filename.empty()) return; if(".ffp"!=filename.substr(filename.size()-4, 4)) filename += ".ffp"; boost::property_tree::ptree parameters; // fiber generation parameters - parameters.put("fiberfox.fibers.realtime", m_FiberGenerationParameters.m_RealTimeFibers); - parameters.put("fiberfox.fibers.showadvanced", m_FiberGenerationParameters.m_AdvancedOptions); - parameters.put("fiberfox.fibers.distribution", m_FiberGenerationParameters.m_Distribution); - parameters.put("fiberfox.fibers.variance", m_FiberGenerationParameters.m_Variance); - parameters.put("fiberfox.fibers.density", m_FiberGenerationParameters.m_FiberDensity); - parameters.put("fiberfox.fibers.spline.sampling", m_FiberGenerationParameters.m_Sampling); - parameters.put("fiberfox.fibers.spline.tension", m_FiberGenerationParameters.m_Tension); - parameters.put("fiberfox.fibers.spline.continuity", m_FiberGenerationParameters.m_Continuity); - parameters.put("fiberfox.fibers.spline.bias", m_FiberGenerationParameters.m_Bias); - parameters.put("fiberfox.fibers.constantradius", m_FiberGenerationParameters.m_ConstantRadius); - parameters.put("fiberfox.fibers.rotation.x", m_FiberGenerationParameters.m_Rotation[0]); - parameters.put("fiberfox.fibers.rotation.y", m_FiberGenerationParameters.m_Rotation[1]); - parameters.put("fiberfox.fibers.rotation.z", m_FiberGenerationParameters.m_Rotation[2]); - parameters.put("fiberfox.fibers.translation.x", m_FiberGenerationParameters.m_Translation[0]); - parameters.put("fiberfox.fibers.translation.y", m_FiberGenerationParameters.m_Translation[1]); - parameters.put("fiberfox.fibers.translation.z", m_FiberGenerationParameters.m_Translation[2]); - parameters.put("fiberfox.fibers.scale.x", m_FiberGenerationParameters.m_Scale[0]); - parameters.put("fiberfox.fibers.scale.y", m_FiberGenerationParameters.m_Scale[1]); - parameters.put("fiberfox.fibers.scale.z", m_FiberGenerationParameters.m_Scale[2]); - parameters.put("fiberfox.fibers.includeFiducials", m_FiberGenerationParameters.m_IncludeFiducials); + parameters.put("fiberfox.fibers.realtime", m_FiberGen.m_RealTimeFibers); + parameters.put("fiberfox.fibers.showadvanced", m_FiberGen.m_AdvancedOptions); + 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_FiberDensity); + 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.constantradius", m_FiberGen.m_ConstantRadius); + 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]); + parameters.put("fiberfox.fibers.includeFiducials", m_FiberGen.m_IncludeFiducials); // image generation parameters - parameters.put("fiberfox.image.basic.size.x", m_ImageRegion.GetSize(0)); - parameters.put("fiberfox.image.basic.size.y", m_ImageRegion.GetSize(1)); - parameters.put("fiberfox.image.basic.size.z", m_ImageRegion.GetSize(2)); - parameters.put("fiberfox.image.basic.spacing.x", m_ImageSpacing[0]); - parameters.put("fiberfox.image.basic.spacing.y", m_ImageSpacing[1]); - parameters.put("fiberfox.image.basic.spacing.z", m_ImageSpacing[2]); - parameters.put("fiberfox.image.basic.origin.x", m_ImageOrigin[0]); - parameters.put("fiberfox.image.basic.origin.y", m_ImageOrigin[1]); - parameters.put("fiberfox.image.basic.origin.z", m_ImageOrigin[2]); - parameters.put("fiberfox.image.basic.direction.1", m_ImageDirection[0][0]); - parameters.put("fiberfox.image.basic.direction.2", m_ImageDirection[0][1]); - parameters.put("fiberfox.image.basic.direction.3", m_ImageDirection[0][2]); - parameters.put("fiberfox.image.basic.direction.4", m_ImageDirection[1][0]); - parameters.put("fiberfox.image.basic.direction.5", m_ImageDirection[1][1]); - parameters.put("fiberfox.image.basic.direction.6", m_ImageDirection[1][2]); - parameters.put("fiberfox.image.basic.direction.7", m_ImageDirection[2][0]); - parameters.put("fiberfox.image.basic.direction.8", m_ImageDirection[2][1]); - parameters.put("fiberfox.image.basic.direction.9", m_ImageDirection[2][2]); - - parameters.put("fiberfox.image.signalScale", m_SignalScale); - parameters.put("fiberfox.image.tEcho", m_tEcho); - parameters.put("fiberfox.image.tLine", m_tLine); - parameters.put("fiberfox.image.tInhom", m_tInhom); - parameters.put("fiberfox.image.bvalue", m_Bvalue); - parameters.put("fiberfox.image.simulatekspace", m_SimulateKspaceAcquisition); - - parameters.put("fiberfox.image.axonRadius", m_AxonRadius); - parameters.put("fiberfox.image.diffusiondirectionmode", m_DiffusionDirectionMode); - parameters.put("fiberfox.image.fiberseparationthreshold", m_FiberSeparationThreshold); - - parameters.put("fiberfox.image.artifacts.addnoise", m_DoAddNoise); - parameters.put("fiberfox.image.artifacts.addghosts", m_DoAddGhosts); - parameters.put("fiberfox.image.artifacts.addaliasing", m_DoAddAliasing); - parameters.put("fiberfox.image.artifacts.addspikes", m_DoAddSpikes); - parameters.put("fiberfox.image.artifacts.addeddycurrents", m_DoAddEddyCurrents); - parameters.put("fiberfox.image.artifacts.spikesnum", m_Spikes); - parameters.put("fiberfox.image.artifacts.spikesscale", m_SpikeAmplitude); - parameters.put("fiberfox.image.artifacts.kspaceLineOffset", m_KspaceLineOffset); - parameters.put("fiberfox.image.artifacts.eddyStrength", m_EddyStrength); - parameters.put("fiberfox.image.artifacts.eddyTau", m_Tau); - parameters.put("fiberfox.image.artifacts.aliasingfactor", m_CroppingFactor); - parameters.put("fiberfox.image.artifacts.addringing", m_DoAddGibbsRinging); - parameters.put("fiberfox.image.doSimulateRelaxation", m_DoSimulateRelaxation); - parameters.put("fiberfox.image.doDisablePartialVolume", m_DoDisablePartialVolume); - parameters.put("fiberfox.image.artifacts.doAddMotion", m_DoAddMotion); - parameters.put("fiberfox.image.artifacts.randomMotion", m_DoRandomizeMotion); - parameters.put("fiberfox.image.artifacts.doAddDistortions", m_DoAddDistortions); - parameters.put("fiberfox.image.artifacts.translation0", m_Translation[0]); - parameters.put("fiberfox.image.artifacts.translation1", m_Translation[1]); - parameters.put("fiberfox.image.artifacts.translation2", m_Translation[2]); - parameters.put("fiberfox.image.artifacts.rotation0", m_Rotation[0]); - parameters.put("fiberfox.image.artifacts.rotation1", m_Rotation[1]); - parameters.put("fiberfox.image.artifacts.rotation2", m_Rotation[2]); - - parameters.put("fiberfox.image.signalmodelstring", m_SignalModelString); - parameters.put("fiberfox.image.artifactmodelstring", m_ArtifactModelString); - parameters.put("fiberfox.image.outpath", m_OutputPath); - parameters.put("fiberfox.image.outputvolumefractions", m_OutputVolumeFractions); - parameters.put("fiberfox.image.showadvanced", m_AdvancedOptions); + 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.1", m_SignalGen.m_ImageDirection[0][0]); + parameters.put("fiberfox.image.basic.direction.2", m_SignalGen.m_ImageDirection[0][1]); + parameters.put("fiberfox.image.basic.direction.3", m_SignalGen.m_ImageDirection[0][2]); + parameters.put("fiberfox.image.basic.direction.4", m_SignalGen.m_ImageDirection[1][0]); + parameters.put("fiberfox.image.basic.direction.5", m_SignalGen.m_ImageDirection[1][1]); + parameters.put("fiberfox.image.basic.direction.6", m_SignalGen.m_ImageDirection[1][2]); + parameters.put("fiberfox.image.basic.direction.7", m_SignalGen.m_ImageDirection[2][0]); + parameters.put("fiberfox.image.basic.direction.8", m_SignalGen.m_ImageDirection[2][1]); + parameters.put("fiberfox.image.basic.direction.9", m_SignalGen.m_ImageDirection[2][2]); + + 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.bvalue", m_SignalGen.m_Bvalue); + parameters.put("fiberfox.image.simulatekspace", m_SignalGen.m_SimulateKspaceAcquisition); + + parameters.put("fiberfox.image.axonRadius", m_SignalGen.m_AxonRadius); + parameters.put("fiberfox.image.diffusiondirectionmode", m_SignalGen.m_DiffusionDirectionMode); + parameters.put("fiberfox.image.fiberseparationthreshold", m_SignalGen.m_FiberSeparationThreshold); + + parameters.put("fiberfox.image.artifacts.addnoise", m_SignalGen.m_DoAddNoise); + parameters.put("fiberfox.image.artifacts.addghosts", m_SignalGen.m_DoAddGhosts); + parameters.put("fiberfox.image.artifacts.addaliasing", m_SignalGen.m_DoAddAliasing); + parameters.put("fiberfox.image.artifacts.addspikes", m_SignalGen.m_DoAddSpikes); + parameters.put("fiberfox.image.artifacts.addeddycurrents", m_SignalGen.m_DoAddEddyCurrents); + 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.addringing", m_SignalGen.m_DoAddGibbsRinging); + parameters.put("fiberfox.image.doSimulateRelaxation", m_SignalGen.m_DoSimulateRelaxation); + parameters.put("fiberfox.image.doDisablePartialVolume", m_SignalGen.m_DoDisablePartialVolume); + 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.doAddDistortions", m_SignalGen.m_DoAddDistortions); + 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.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.image.outputvolumefractions", m_Misc.m_OutputVolumeFractions); + parameters.put("fiberfox.image.showadvanced", m_Misc.m_AdvancedOptions); parameters.put("fiberfox.image.basic.numgradients", m_NumGradients); if (m_NoiseModel!=NULL) { parameters.put("fiberfox.image.artifacts.noisevariance", m_NoiseModel->GetNoiseVariance()); if (dynamic_cast*>(m_NoiseModel)) parameters.put("fiberfox.image.artifacts.noisetype", "rice"); else if (dynamic_cast*>(m_NoiseModel)) parameters.put("fiberfox.image.artifacts.noisetype", "chisquare"); } for (int i=0; i* signalModel = NULL; if (i(i)+".type", "fiber"); } else { signalModel = m_NonFiberModelList.at(i-m_FiberModelList.size()); parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".type", "non-fiber"); } if (dynamic_cast*>(signalModel)) { mitk::StickModel* model = dynamic_cast*>(signalModel); parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".stick.d", model->GetDiffusivity()); parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".stick.t2", model->GetT2()); } else if (dynamic_cast*>(signalModel)) { mitk::TensorModel* model = dynamic_cast*>(signalModel); parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".tensor.d1", model->GetDiffusivity1()); parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".tensor.d2", model->GetDiffusivity2()); parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".tensor.d3", model->GetDiffusivity3()); parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".tensor.t2", model->GetT2()); } else if (dynamic_cast*>(signalModel)) { mitk::RawShModel* model = dynamic_cast*>(signalModel); parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".prototype.minFA", model->GetFaRange().first); parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".prototype.maxFA", model->GetFaRange().second); parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".prototype.minADC", model->GetAdcRange().first); parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".prototype.maxADC", model->GetAdcRange().second); parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".prototype.numSamples", model->GetMaxNumKernels()); } else if (dynamic_cast*>(signalModel)) { mitk::BallModel* model = dynamic_cast*>(signalModel); parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".ball.d", model->GetDiffusivity()); parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".ball.t2", model->GetT2()); } else if (dynamic_cast*>(signalModel)) { mitk::AstroStickModel* model = dynamic_cast*>(signalModel); parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".astrosticks.d", model->GetDiffusivity()); parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".astrosticks.t2", model->GetT2()); parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".astrosticks.randomize", model->GetRandomizeSticks()); } else if (dynamic_cast*>(signalModel)) { mitk::DotModel* model = dynamic_cast*>(signalModel); parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".dot.t2", model->GetT2()); } parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".ID", signalModel->m_CompartmentId); } boost::property_tree::xml_writer_settings writerSettings(' ', 2); boost::property_tree::xml_parser::write_xml(filename, parameters, std::locale(), writerSettings); } template< class ScalarType > void mitk::FiberfoxParameters< ScalarType >::LoadParameters(string filename) { if(filename.empty()) return; boost::property_tree::ptree parameterTree; boost::property_tree::xml_parser::read_xml(filename, parameterTree); m_FiberModelList.clear(); m_NonFiberModelList.clear(); if (m_NoiseModel!=NULL) delete m_NoiseModel; BOOST_FOREACH( boost::property_tree::ptree::value_type const& v1, parameterTree.get_child("fiberfox") ) { if( v1.first == "fibers" ) { - m_FiberGenerationParameters.m_RealTimeFibers = v1.second.get("realtime", m_FiberGenerationParameters.m_RealTimeFibers); - m_FiberGenerationParameters.m_AdvancedOptions = v1.second.get("showadvanced", m_FiberGenerationParameters.m_AdvancedOptions); - m_FiberGenerationParameters.m_Distribution = v1.second.get("distribution", m_FiberGenerationParameters.m_Distribution); - m_FiberGenerationParameters.m_Variance = v1.second.get("variance", m_FiberGenerationParameters.m_Variance); - m_FiberGenerationParameters.m_FiberDensity = v1.second.get("density", m_FiberGenerationParameters.m_FiberDensity); - m_FiberGenerationParameters.m_Sampling = v1.second.get("spline.sampling", m_FiberGenerationParameters.m_Sampling); - m_FiberGenerationParameters.m_Tension = v1.second.get("spline.tension", m_FiberGenerationParameters.m_Tension); - m_FiberGenerationParameters.m_Continuity = v1.second.get("spline.continuity", m_FiberGenerationParameters.m_Continuity); - m_FiberGenerationParameters.m_Bias = v1.second.get("spline.bias", m_FiberGenerationParameters.m_Bias); - m_FiberGenerationParameters.m_ConstantRadius = v1.second.get("constantradius", m_FiberGenerationParameters.m_ConstantRadius); - m_FiberGenerationParameters.m_Rotation[0] = v1.second.get("rotation.x", m_FiberGenerationParameters.m_Rotation[0]); - m_FiberGenerationParameters.m_Rotation[1] = v1.second.get("rotation.y", m_FiberGenerationParameters.m_Rotation[1]); - m_FiberGenerationParameters.m_Rotation[2] = v1.second.get("rotation.z", m_FiberGenerationParameters.m_Rotation[2]); - m_FiberGenerationParameters.m_Translation[0] = v1.second.get("translation.x", m_FiberGenerationParameters.m_Translation[0]); - m_FiberGenerationParameters.m_Translation[1] = v1.second.get("translation.y", m_FiberGenerationParameters.m_Translation[1]); - m_FiberGenerationParameters.m_Translation[2] = v1.second.get("translation.z", m_FiberGenerationParameters.m_Translation[2]); - m_FiberGenerationParameters.m_Scale[0] = v1.second.get("scale.x", m_FiberGenerationParameters.m_Scale[0]); - m_FiberGenerationParameters.m_Scale[1] = v1.second.get("scale.y", m_FiberGenerationParameters.m_Scale[1]); - m_FiberGenerationParameters.m_Scale[2] = v1.second.get("scale.z", m_FiberGenerationParameters.m_Scale[2]); - m_FiberGenerationParameters.m_IncludeFiducials = v1.second.get("includeFiducials", m_FiberGenerationParameters.m_IncludeFiducials); + m_FiberGen.m_RealTimeFibers = v1.second.get("realtime", m_FiberGen.m_RealTimeFibers); + m_FiberGen.m_AdvancedOptions = v1.second.get("showadvanced", m_FiberGen.m_AdvancedOptions); + m_FiberGen.m_Distribution = v1.second.get("distribution", m_FiberGen.m_Distribution); + m_FiberGen.m_Variance = v1.second.get("variance", m_FiberGen.m_Variance); + m_FiberGen.m_FiberDensity = v1.second.get("density", m_FiberGen.m_FiberDensity); + m_FiberGen.m_Sampling = v1.second.get("spline.sampling", m_FiberGen.m_Sampling); + m_FiberGen.m_Tension = v1.second.get("spline.tension", m_FiberGen.m_Tension); + m_FiberGen.m_Continuity = v1.second.get("spline.continuity", m_FiberGen.m_Continuity); + m_FiberGen.m_Bias = v1.second.get("spline.bias", m_FiberGen.m_Bias); + m_FiberGen.m_ConstantRadius = v1.second.get("constantradius", m_FiberGen.m_ConstantRadius); + m_FiberGen.m_Rotation[0] = v1.second.get("rotation.x", m_FiberGen.m_Rotation[0]); + m_FiberGen.m_Rotation[1] = v1.second.get("rotation.y", m_FiberGen.m_Rotation[1]); + m_FiberGen.m_Rotation[2] = v1.second.get("rotation.z", m_FiberGen.m_Rotation[2]); + m_FiberGen.m_Translation[0] = v1.second.get("translation.x", m_FiberGen.m_Translation[0]); + m_FiberGen.m_Translation[1] = v1.second.get("translation.y", m_FiberGen.m_Translation[1]); + m_FiberGen.m_Translation[2] = v1.second.get("translation.z", m_FiberGen.m_Translation[2]); + m_FiberGen.m_Scale[0] = v1.second.get("scale.x", m_FiberGen.m_Scale[0]); + m_FiberGen.m_Scale[1] = v1.second.get("scale.y", m_FiberGen.m_Scale[1]); + m_FiberGen.m_Scale[2] = v1.second.get("scale.z", m_FiberGen.m_Scale[2]); + m_FiberGen.m_IncludeFiducials = v1.second.get("includeFiducials", m_FiberGen.m_IncludeFiducials); } else if ( v1.first == "image" ) { - m_ImageRegion.SetSize(0, v1.second.get("basic.size.x", m_ImageRegion.GetSize(0))); - m_ImageRegion.SetSize(1, v1.second.get("basic.size.y", m_ImageRegion.GetSize(1))); - m_ImageRegion.SetSize(2, v1.second.get("basic.size.z", m_ImageRegion.GetSize(2))); - m_ImageSpacing[0] = v1.second.get("basic.spacing.x", m_ImageSpacing[0]); - m_ImageSpacing[1] = v1.second.get("basic.spacing.y", m_ImageSpacing[1]); - m_ImageSpacing[2] = v1.second.get("basic.spacing.z", m_ImageSpacing[2]); - m_ImageOrigin[0] = v1.second.get("basic.origin.x", m_ImageOrigin[0]); - m_ImageOrigin[1] = v1.second.get("basic.origin.y", m_ImageOrigin[1]); - m_ImageOrigin[2] = v1.second.get("basic.origin.z", m_ImageOrigin[2]); - m_ImageDirection[0][0] = v1.second.get("basic.direction.1", m_ImageDirection[0][0]); - m_ImageDirection[0][1] = v1.second.get("basic.direction.2", m_ImageDirection[0][1]); - m_ImageDirection[0][2] = v1.second.get("basic.direction.3", m_ImageDirection[0][2]); - m_ImageDirection[1][0] = v1.second.get("basic.direction.4", m_ImageDirection[1][0]); - m_ImageDirection[1][1] = v1.second.get("basic.direction.5", m_ImageDirection[1][1]); - m_ImageDirection[1][2] = v1.second.get("basic.direction.6", m_ImageDirection[1][2]); - m_ImageDirection[2][0] = v1.second.get("basic.direction.7", m_ImageDirection[2][0]); - m_ImageDirection[2][1] = v1.second.get("basic.direction.8", m_ImageDirection[2][1]); - m_ImageDirection[2][2] = v1.second.get("basic.direction.9", m_ImageDirection[2][2]); - - m_SignalScale = v1.second.get("signalScale", m_SignalScale); - m_tEcho = v1.second.get("tEcho", m_tEcho); - m_tLine = v1.second.get("tLine", m_tLine); - m_tInhom = v1.second.get("tInhom", m_tInhom); - m_Bvalue = v1.second.get("bvalue", m_Bvalue); - m_SimulateKspaceAcquisition = v1.second.get("simulatekspace", m_SimulateKspaceAcquisition); - - m_AxonRadius = v1.second.get("axonRadius", m_AxonRadius); + m_SignalGen.m_ImageRegion.SetSize(0, v1.second.get("basic.size.x",m_SignalGen.m_ImageRegion.GetSize(0))); + m_SignalGen.m_ImageRegion.SetSize(1, v1.second.get("basic.size.y",m_SignalGen.m_ImageRegion.GetSize(1))); + m_SignalGen.m_ImageRegion.SetSize(2, v1.second.get("basic.size.z",m_SignalGen.m_ImageRegion.GetSize(2))); + m_SignalGen.m_ImageSpacing[0] = v1.second.get("basic.spacing.x",m_SignalGen.m_ImageSpacing[0]); + m_SignalGen.m_ImageSpacing[1] = v1.second.get("basic.spacing.y",m_SignalGen.m_ImageSpacing[1]); + m_SignalGen.m_ImageSpacing[2] = v1.second.get("basic.spacing.z",m_SignalGen.m_ImageSpacing[2]); + m_SignalGen.m_ImageOrigin[0] = v1.second.get("basic.origin.x",m_SignalGen.m_ImageOrigin[0]); + m_SignalGen.m_ImageOrigin[1] = v1.second.get("basic.origin.y",m_SignalGen.m_ImageOrigin[1]); + m_SignalGen.m_ImageOrigin[2] = v1.second.get("basic.origin.z",m_SignalGen.m_ImageOrigin[2]); + m_SignalGen.m_ImageDirection[0][0] = v1.second.get("basic.direction.1",m_SignalGen.m_ImageDirection[0][0]); + m_SignalGen.m_ImageDirection[0][1] = v1.second.get("basic.direction.2",m_SignalGen.m_ImageDirection[0][1]); + m_SignalGen.m_ImageDirection[0][2] = v1.second.get("basic.direction.3",m_SignalGen.m_ImageDirection[0][2]); + m_SignalGen.m_ImageDirection[1][0] = v1.second.get("basic.direction.4",m_SignalGen.m_ImageDirection[1][0]); + m_SignalGen.m_ImageDirection[1][1] = v1.second.get("basic.direction.5",m_SignalGen.m_ImageDirection[1][1]); + m_SignalGen.m_ImageDirection[1][2] = v1.second.get("basic.direction.6",m_SignalGen.m_ImageDirection[1][2]); + m_SignalGen.m_ImageDirection[2][0] = v1.second.get("basic.direction.7",m_SignalGen.m_ImageDirection[2][0]); + m_SignalGen.m_ImageDirection[2][1] = v1.second.get("basic.direction.8",m_SignalGen.m_ImageDirection[2][1]); + m_SignalGen.m_ImageDirection[2][2] = v1.second.get("basic.direction.9",m_SignalGen.m_ImageDirection[2][2]); + + m_SignalGen.m_SignalScale = v1.second.get("signalScale", m_SignalGen.m_SignalScale); + m_SignalGen.m_tEcho = v1.second.get("tEcho", m_SignalGen.m_tEcho); + m_SignalGen.m_tLine = v1.second.get("tLine", m_SignalGen.m_tLine); + m_SignalGen.m_tInhom = v1.second.get("tInhom", m_SignalGen.m_tInhom); + m_SignalGen.m_Bvalue = v1.second.get("bvalue", m_SignalGen.m_Bvalue); + m_SignalGen.m_SimulateKspaceAcquisition = v1.second.get("simulatekspace", m_SignalGen.m_SimulateKspaceAcquisition); + + m_SignalGen.m_AxonRadius = v1.second.get("axonRadius", m_SignalGen.m_AxonRadius); int mode = v1.second.get("diffusiondirectionmode", 0); switch (mode) { case 0: - m_DiffusionDirectionMode = FIBER_TANGENT_DIRECTIONS; + m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::FIBER_TANGENT_DIRECTIONS; break; case 1: - m_DiffusionDirectionMode = MAIN_FIBER_DIRECTIONS; + m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::MAIN_FIBER_DIRECTIONS; break; case 2: - m_DiffusionDirectionMode = RANDOM_DIRECTIONS; + m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::RANDOM_DIRECTIONS; break; default: - m_DiffusionDirectionMode = FIBER_TANGENT_DIRECTIONS; + m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::FIBER_TANGENT_DIRECTIONS; } - m_FiberSeparationThreshold = v1.second.get("fiberseparationthreshold", m_FiberSeparationThreshold); - - m_DoAddNoise = v1.second.get("artifacts.addnoise", m_DoAddNoise); - m_DoAddGhosts = v1.second.get("artifacts.addghosts", m_DoAddGhosts); - m_DoAddAliasing = v1.second.get("artifacts.addaliasing", m_DoAddAliasing); - m_DoAddSpikes = v1.second.get("artifacts.addspikes", m_DoAddSpikes); - m_DoAddEddyCurrents = v1.second.get("artifacts.addeddycurrents", m_DoAddEddyCurrents); - m_Spikes = v1.second.get("artifacts.spikesnum", m_Spikes); - m_SpikeAmplitude = v1.second.get("artifacts.spikesscale", m_SpikeAmplitude); - m_KspaceLineOffset = v1.second.get("artifacts.kspaceLineOffset", m_KspaceLineOffset); - m_EddyStrength = v1.second.get("artifacts.eddyStrength", m_EddyStrength); - m_Tau = v1.second.get("artifacts.eddyTau", m_Tau); - m_CroppingFactor = v1.second.get("artifacts.aliasingfactor", m_CroppingFactor); - m_DoAddGibbsRinging = v1.second.get("artifacts.addringing", m_DoAddGibbsRinging); - m_DoSimulateRelaxation = v1.second.get("doSimulateRelaxation", m_DoSimulateRelaxation); - m_DoDisablePartialVolume = v1.second.get("doDisablePartialVolume", m_DoDisablePartialVolume); - m_DoAddMotion = v1.second.get("artifacts.doAddMotion", m_DoAddMotion); - m_DoRandomizeMotion = v1.second.get("artifacts.randomMotion", m_DoRandomizeMotion); - m_DoAddDistortions = v1.second.get("artifacts.doAddDistortions", m_DoAddDistortions); - m_Translation[0] = v1.second.get("artifacts.translation0", m_Translation[0]); - m_Translation[1] = v1.second.get("artifacts.translation1", m_Translation[1]); - m_Translation[2] = v1.second.get("artifacts.translation2", m_Translation[2]); - m_Rotation[0] = v1.second.get("artifacts.rotation0", m_Rotation[0]); - m_Rotation[1] = v1.second.get("artifacts.rotation1", m_Rotation[1]); - m_Rotation[2] = v1.second.get("artifacts.rotation2", m_Rotation[2]); - - m_SignalModelString = v1.second.get("signalmodelstring", m_SignalModelString); - m_ArtifactModelString = v1.second.get("artifactmodelstring", m_ArtifactModelString); - m_OutputPath = v1.second.get("outpath", m_OutputPath); - m_OutputVolumeFractions = v1.second.get("outputvolumefractions", m_OutputVolumeFractions); - m_AdvancedOptions = v1.second.get("showadvanced", m_AdvancedOptions); + m_SignalGen.m_FiberSeparationThreshold = v1.second.get("fiberseparationthreshold", m_SignalGen.m_FiberSeparationThreshold); + + m_SignalGen.m_DoAddNoise = v1.second.get("artifacts.addnoise", m_SignalGen.m_DoAddNoise); + m_SignalGen.m_DoAddGhosts = v1.second.get("artifacts.addghosts", m_SignalGen.m_DoAddGhosts); + m_SignalGen.m_DoAddAliasing = v1.second.get("artifacts.addaliasing", m_SignalGen.m_DoAddAliasing); + m_SignalGen.m_DoAddSpikes = v1.second.get("artifacts.addspikes", m_SignalGen.m_DoAddSpikes); + m_SignalGen.m_DoAddEddyCurrents = v1.second.get("artifacts.addeddycurrents", m_SignalGen.m_DoAddEddyCurrents); + m_SignalGen.m_Spikes = v1.second.get("artifacts.spikesnum", m_SignalGen.m_Spikes); + m_SignalGen.m_SpikeAmplitude = v1.second.get("artifacts.spikesscale", m_SignalGen.m_SpikeAmplitude); + m_SignalGen.m_KspaceLineOffset = v1.second.get("artifacts.kspaceLineOffset", m_SignalGen.m_KspaceLineOffset); + m_SignalGen.m_EddyStrength = v1.second.get("artifacts.eddyStrength", m_SignalGen.m_EddyStrength); + m_SignalGen.m_Tau = v1.second.get("artifacts.eddyTau", m_SignalGen.m_Tau); + m_SignalGen.m_CroppingFactor = v1.second.get("artifacts.aliasingfactor", m_SignalGen.m_CroppingFactor); + m_SignalGen.m_DoAddGibbsRinging = v1.second.get("artifacts.addringing", m_SignalGen.m_DoAddGibbsRinging); + m_SignalGen.m_DoSimulateRelaxation = v1.second.get("doSimulateRelaxation", m_SignalGen.m_DoSimulateRelaxation); + m_SignalGen.m_DoDisablePartialVolume = v1.second.get("doDisablePartialVolume", m_SignalGen.m_DoDisablePartialVolume); + m_SignalGen.m_DoAddMotion = v1.second.get("artifacts.doAddMotion", m_SignalGen.m_DoAddMotion); + m_SignalGen.m_DoRandomizeMotion = v1.second.get("artifacts.randomMotion", m_SignalGen.m_DoRandomizeMotion); + m_SignalGen.m_DoAddDistortions = v1.second.get("artifacts.doAddDistortions", m_SignalGen.m_DoAddDistortions); + m_SignalGen.m_Translation[0] = v1.second.get("artifacts.translation0", m_SignalGen.m_Translation[0]); + m_SignalGen.m_Translation[1] = v1.second.get("artifacts.translation1", m_SignalGen.m_Translation[1]); + m_SignalGen.m_Translation[2] = v1.second.get("artifacts.translation2", m_SignalGen.m_Translation[2]); + m_SignalGen.m_Rotation[0] = v1.second.get("artifacts.rotation0", m_SignalGen.m_Rotation[0]); + m_SignalGen.m_Rotation[1] = v1.second.get("artifacts.rotation1", m_SignalGen.m_Rotation[1]); + m_SignalGen.m_Rotation[2] = v1.second.get("artifacts.rotation2", m_SignalGen.m_Rotation[2]); + + m_Misc.m_SignalModelString = v1.second.get("signalmodelstring", m_Misc.m_SignalModelString); + m_Misc.m_ArtifactModelString = v1.second.get("artifactmodelstring", m_Misc.m_ArtifactModelString); + m_Misc.m_OutputPath = v1.second.get("outpath", m_Misc.m_OutputPath); + m_Misc.m_OutputVolumeFractions = v1.second.get("outputvolumefractions", m_Misc.m_OutputVolumeFractions); + m_Misc.m_AdvancedOptions = v1.second.get("showadvanced", m_Misc.m_AdvancedOptions); m_NumGradients = v1.second.get("numgradients", m_NumGradients); GenerateGradientHalfShell(); try { if (v1.second.get("artifacts.noisetype")=="rice") { m_NoiseModel = new mitk::RicianNoiseModel(); m_NoiseModel->SetNoiseVariance(v1.second.get("artifacts.noisevariance")); } } catch(...) {} try { if (v1.second.get("artifacts.noisetype")=="chisquare") { m_NoiseModel = new mitk::ChiSquareNoiseModel(); m_NoiseModel->SetNoiseVariance(v1.second.get("artifacts.noisevariance")); } } catch(...){ } BOOST_FOREACH( boost::property_tree::ptree::value_type const& v2, v1.second.get_child("compartments") ) { try // stick model { mitk::StickModel* model = new mitk::StickModel(); model->SetDiffusivity(v2.second.get("stick.d")); model->SetT2(v2.second.get("stick.t2")); model->m_CompartmentId = v2.second.get("ID"); if (v2.second.get("type")=="fiber") m_FiberModelList.push_back(model); else if (v2.second.get("type")=="non-fiber") m_NonFiberModelList.push_back(model); } catch (...) { } try // tensor model { mitk::TensorModel* model = new mitk::TensorModel(); model->SetDiffusivity1(v2.second.get("tensor.d1")); model->SetDiffusivity2(v2.second.get("tensor.d2")); model->SetDiffusivity3(v2.second.get("tensor.d3")); model->SetT2(v2.second.get("tensor.t2")); model->m_CompartmentId = v2.second.get("ID"); if (v2.second.get("type")=="fiber") m_FiberModelList.push_back(model); else if (v2.second.get("type")=="non-fiber") m_NonFiberModelList.push_back(model); } catch (...) { } try // ball model { mitk::BallModel* model = new mitk::BallModel(); model->SetDiffusivity(v2.second.get("ball.d")); model->SetT2(v2.second.get("ball.t2")); model->m_CompartmentId = v2.second.get("ID"); if (v2.second.get("type")=="fiber") m_FiberModelList.push_back(model); else if (v2.second.get("type")=="non-fiber") m_NonFiberModelList.push_back(model); } catch (...) { } try // astrosticks model { mitk::AstroStickModel* model = new AstroStickModel(); model->SetDiffusivity(v2.second.get("astrosticks.d")); model->SetT2(v2.second.get("astrosticks.t2")); model->SetRandomizeSticks(v2.second.get("astrosticks.randomize")); model->m_CompartmentId = v2.second.get("ID"); if (v2.second.get("type")=="fiber") m_FiberModelList.push_back(model); else if (v2.second.get("type")=="non-fiber") m_NonFiberModelList.push_back(model); } catch (...) { } try // dot model { mitk::DotModel* model = new mitk::DotModel(); model->SetT2(v2.second.get("dot.t2")); model->m_CompartmentId = v2.second.get("ID"); if (v2.second.get("type")=="fiber") m_FiberModelList.push_back(model); else if (v2.second.get("type")=="non-fiber") m_NonFiberModelList.push_back(model); } catch (...) { } try // prototype model { mitk::RawShModel* model = new mitk::RawShModel(); model->SetMaxNumKernels(v2.second.get("prototype.numSamples")); model->SetFaRange(v2.second.get("prototype.minFA"), v2.second.get("prototype.maxFA")); model->SetAdcRange(v2.second.get("prototype.minADC"), v2.second.get("prototype.maxADC")); model->m_CompartmentId = v2.second.get("ID"); if (v2.second.get("type")=="fiber") m_FiberModelList.push_back(model); else if (v2.second.get("type")=="non-fiber") m_NonFiberModelList.push_back(model); } catch (...) { } } } } } template< class ScalarType > void mitk::FiberfoxParameters< ScalarType >::PrintSelf() { MITK_INFO << "Not implemented :("; } diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h index 2359106647..28f4aa0ce3 100644 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h @@ -1,219 +1,206 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _MITK_FiberfoxParameters_H #define _MITK_FiberfoxParameters_H #include #include #include #include #include #include #include #include using namespace std; namespace mitk { -/** - * \brief Datastructure to manage the Fiberfox signal generation parameters. - * - */ +/** Fiber generation */ +struct FiberGenerationParameters +{ + bool m_RealTimeFibers; + bool m_AdvancedOptions; + int m_Distribution; + double m_Variance; + unsigned int m_FiberDensity; + double m_Sampling; + double m_Tension; + double m_Continuity; + double m_Bias; + bool m_ConstantRadius; + mitk::Vector3D m_Rotation; + mitk::Vector3D m_Translation; + mitk::Vector3D m_Scale; + bool m_IncludeFiducials; +}; -template< class ScalarType = double > -class FiberfoxParameters +/** Signal generation */ +struct SignalGenerationParameters { -public: + typedef itk::Image ItkDoubleImgType; + typedef itk::Image ItkUcharImgType; enum DiffusionDirectionMode { FIBER_TANGENT_DIRECTIONS, MAIN_FIBER_DIRECTIONS, RANDOM_DIRECTIONS }; - typedef itk::Image ItkDoubleImgType; - typedef itk::Image ItkUcharImgType; - typedef std::vector< DiffusionSignalModel* > DiffusionModelListType; - typedef DiffusionSignalModel::GradientListType GradientListType; - typedef DiffusionSignalModel::GradientType GradientType; - typedef DiffusionNoiseModel NoiseModelType; - typedef DiffusionSignalModel* DiffusionModelType; - - FiberfoxParameters(); - ~FiberfoxParameters(); - - /** Get same parameter object with different template parameter */ - template< class OutType > FiberfoxParameters< OutType > CopyParameters() - { - FiberfoxParameters< OutType > out; - -// out.m_FiberGenerationParameters = m_FiberGenerationParameters; - - out.m_ImageRegion = m_ImageRegion; - out.m_ImageSpacing = m_ImageSpacing; - out.m_ImageOrigin = m_ImageOrigin; - out.m_ImageDirection = m_ImageDirection; - out.SetNumWeightedGradients(m_NumGradients); - out.m_Bvalue = m_Bvalue; - out.m_SignalScale = m_SignalScale; - out.m_tEcho = m_tEcho; - out.m_tLine = m_tLine; - out.m_tInhom = m_tInhom; - out.m_AxonRadius = m_AxonRadius; - out.m_KspaceLineOffset = m_KspaceLineOffset; - out.m_DoAddGibbsRinging = m_DoAddGibbsRinging; - out.m_EddyStrength = m_EddyStrength; - out.m_Spikes = m_Spikes; - out.m_SpikeAmplitude = m_SpikeAmplitude; - out.m_CroppingFactor = m_CroppingFactor; - out.m_DoSimulateRelaxation = m_DoSimulateRelaxation; - out.m_DoDisablePartialVolume = m_DoDisablePartialVolume; - out.m_DoAddMotion = m_DoAddMotion; - out.m_DoRandomizeMotion = m_DoRandomizeMotion; - out.m_Translation = m_Translation; - out.m_Rotation = m_Rotation; - if (m_NoiseModel!=NULL) - { - if (dynamic_cast*>(m_NoiseModel)) - out.m_NoiseModel = new mitk::RicianNoiseModel(); - else if (dynamic_cast*>(m_NoiseModel)) - out.m_NoiseModel = new mitk::ChiSquareNoiseModel(); - out.m_NoiseModel->SetNoiseVariance(m_NoiseModel->GetNoiseVariance()); - } - out.m_FrequencyMap = m_FrequencyMap; - out.m_MaskImage = m_MaskImage; - out.m_ResultNode = m_ResultNode; - out.m_ParentNode = m_ParentNode; - out.m_SignalModelString = m_SignalModelString; - out.m_ArtifactModelString = m_ArtifactModelString; - out.m_OutputPath = m_OutputPath; - // out.m_DiffusionDirectionMode = m_DiffusionDirectionMode; - // out.m_SignalGenerationMode = m_SignalGenerationMode; - out.m_SimulateKspaceAcquisition = m_SimulateKspaceAcquisition; - - // TODO: copy constructor für singalmodelle und rauschen - - return out; - } - - /** Fiber generation */ - struct FiberGenerationParameters - { - bool m_RealTimeFibers; - bool m_AdvancedOptions; - int m_Distribution; - double m_Variance; - unsigned int m_FiberDensity; - double m_Sampling; - double m_Tension; - double m_Continuity; - double m_Bias; - bool m_ConstantRadius; - mitk::Vector3D m_Rotation; - mitk::Vector3D m_Translation; - mitk::Vector3D m_Scale; - bool m_IncludeFiducials; - } m_FiberGenerationParameters; - - /** Output image specifications */ + /** input/output image specifications */ itk::ImageRegion<3> m_ImageRegion; ///< Image size. itk::Vector m_ImageSpacing; ///< Image voxel size. itk::Point m_ImageOrigin; ///< Image origin. itk::Matrix m_ImageDirection; ///< Image rotation matrix. /** Other acquisitions parameters */ double m_SignalScale; ///< Scaling factor for output signal (before noise is added). double m_tEcho; ///< Echo time TE. double m_tLine; ///< k-space line readout time. double m_tInhom; ///< T2' double m_Bvalue; ///< Acquisition b-value bool m_SimulateKspaceAcquisition;///< - /** Signal generation */ - DiffusionModelListType m_FiberModelList; ///< Intra- and inter-axonal compartments. - DiffusionModelListType m_NonFiberModelList; ///< Extra-axonal compartments. double m_AxonRadius; ///< Determines compartment volume fractions (0 == automatic axon radius estimation) DiffusionDirectionMode m_DiffusionDirectionMode; ///< Determines how the main diffusion direction of the signal models is selected double m_FiberSeparationThreshold; ///< Used for random and and mein fiber deriction DiffusionDirectionMode /** Artifacts */ unsigned int m_Spikes; ///< Number of spikes randomly appearing in the image double m_SpikeAmplitude; ///< amplitude of spikes relative to the largest signal intensity (magnitude of complex) double m_KspaceLineOffset; ///< Causes N/2 ghosts. Larger offset means stronger ghost. double m_EddyStrength; ///< Strength of eddy current induced gradients in mT/m. double m_Tau; ///< Eddy current decay constant (in ms) double m_CroppingFactor; ///< FOV size in y-direction is multiplied by this factor. Causes aliasing artifacts. + ItkDoubleImgType::Pointer m_FrequencyMap; ///< If != NULL, distortions are added to the image using this frequency map. + ItkUcharImgType::Pointer m_MaskImage; ///< Signal is only genrated inside of the mask image. + bool m_DoAddNoise; bool m_DoAddGhosts; bool m_DoAddAliasing; bool m_DoAddSpikes; bool m_DoAddEddyCurrents; + bool m_DoAddGibbsRinging; ///< Add Gibbs ringing artifact bool m_DoSimulateRelaxation; ///< Add T2 relaxation effects bool m_DoAddDistortions; ///< Add magnetic field distortions bool m_DoDisablePartialVolume; ///< Disable partial volume effects. Each voxel is either all fiber or all non-fiber. bool m_DoAddMotion; ///< Enable motion artifacts. bool m_DoRandomizeMotion; ///< Toggles between random and linear motion. itk::Vector m_Translation; ///< Maximum translational motion. itk::Vector m_Rotation; ///< Maximum rotational motion. - NoiseModelType* m_NoiseModel; ///< If != NULL, noise is added to the image. - ItkDoubleImgType::Pointer m_FrequencyMap; ///< If != NULL, distortions are added to the image using this frequency map. - ItkUcharImgType::Pointer m_MaskImage; ///< Signal is only genrated inside of the mask image. +}; + +/** GUI persistence, input, output, ... */ +struct MiscFiberfoxParameters +{ + mitk::DataNode::Pointer m_ResultNode; ///< Stores resulting image. + mitk::DataNode::Pointer m_ParentNode; ///< Parent node of result node. + string m_SignalModelString; ///< Appendet to the name of the result node + string m_ArtifactModelString; ///< Appendet to the name of the result node + string m_OutputPath; ///< Image is automatically saved to the specified folder after simulation is finished. + bool m_OutputVolumeFractions; + bool m_AdvancedOptions; +}; + +/** + * \brief Datastructure to manage the Fiberfox signal generation parameters. + * + */ +template< class ScalarType = double > +class FiberfoxParameters +{ +public: + + typedef itk::Image ItkDoubleImgType; + typedef itk::Image ItkUcharImgType; + typedef DiffusionSignalModel DiffusionModelType; + typedef std::vector< DiffusionModelType* > DiffusionModelListType; + typedef typename DiffusionModelType::GradientListType GradientListType; + typedef typename DiffusionModelType::GradientType GradientType; + typedef DiffusionNoiseModel NoiseModelType; + + FiberfoxParameters(); + ~FiberfoxParameters(); + + /** Get same parameter object with different template parameter */ + template< class OutType > FiberfoxParameters< OutType > CopyParameters() + { + FiberfoxParameters< OutType > out; - /** Output parameters (only relevant in GUI application) */ - mitk::DataNode::Pointer m_ResultNode; ///< Stores resulting image. - mitk::DataNode::Pointer m_ParentNode; ///< Parent node of result node. - string m_SignalModelString; ///< Appendet to the name of the result node - string m_ArtifactModelString; ///< Appendet to the name of the result node - string m_OutputPath; ///< Image is automatically saved to the specified folder after simulation is finished. - bool m_OutputVolumeFractions; - bool m_AdvancedOptions; + out.m_FiberGen = m_FiberGen; + out.m_SignalGen = m_SignalGen; + out.m_Misc = m_Misc; + + out.SetNumWeightedGradients(m_NumGradients); + + if (m_NoiseModel!=NULL) + { + if (dynamic_cast*>(m_NoiseModel)) + out.m_NoiseModel = new mitk::RicianNoiseModel(); + else if (dynamic_cast*>(m_NoiseModel)) + out.m_NoiseModel = new mitk::ChiSquareNoiseModel(); + out.m_NoiseModel->SetNoiseVariance(m_NoiseModel->GetNoiseVariance()); + } + + // TODO: copy constructor für singalmodelle, gradienten + + return out; + } + + FiberGenerationParameters m_FiberGen; + SignalGenerationParameters m_SignalGen; + MiscFiberfoxParameters m_Misc; + + /** Signal generation */ + DiffusionModelListType m_FiberModelList; ///< Intra- and inter-axonal compartments. + DiffusionModelListType m_NonFiberModelList; ///< Extra-axonal compartments. + NoiseModelType* m_NoiseModel; ///< If != NULL, noise is added to the image. void PrintSelf(); ///< Print parameters to stdout. void SaveParameters(string filename); ///< Save image generation parameters to .ffp file. void LoadParameters(string filename); ///< Load image generation parameters from .ffp file. void GenerateGradientHalfShell(); ///< Generates half shell of gradient directions (with m_NumGradients non-zero directions) std::vector< int > GetBaselineIndices(); unsigned int GetFirstBaselineIndex(); bool IsBaselineIndex(unsigned int idx); unsigned int GetNumWeightedVolumes(); unsigned int GetNumBaselineVolumes(); unsigned int GetNumVolumes(); GradientListType GetGradientDirections(); GradientType GetGradientDirection(unsigned int i); void SetNumWeightedGradients(int numGradients); ///< Automaticall calls GenerateGradientHalfShell() afterwards. void SetGradienDirections(GradientListType gradientList); void SetGradienDirections(mitk::DiffusionImage::GradientDirectionContainerType::Pointer gradientList); protected: unsigned int m_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. }; } #include "mitkFiberfoxParameters.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxAddArtifactsToDwiTest.cpp b/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxAddArtifactsToDwiTest.cpp index bcb759106d..e170a82360 100644 --- a/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxAddArtifactsToDwiTest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxAddArtifactsToDwiTest.cpp @@ -1,181 +1,181 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include /**Documentation * Test the Fiberfox simulation functions (diffusion weighted image -> diffusion weighted image) */ class mitkFiberfoxAddArtifactsToDwiTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkFiberfoxAddArtifactsToDwiTestSuite); MITK_TEST(Spikes); MITK_TEST(GibbsRinging); MITK_TEST(Ghost); MITK_TEST(Aliasing); MITK_TEST(Eddy); MITK_TEST(RicianNoise); MITK_TEST(ChiSquareNoise); MITK_TEST(Distortions); CPPUNIT_TEST_SUITE_END(); private: mitk::DiffusionImage::Pointer m_InputDwi; FiberfoxParameters m_Parameters; public: void setUp() { // reference files m_InputDwi = dynamic_cast*>(mitk::IOUtil::LoadDataNode(GetTestDataFilePath("DiffusionImaging/Fiberfox/StickBall_RELAX.dwi"))->GetData()); // parameter setup m_Parameters = FiberfoxParameters(); - m_Parameters.m_ImageRegion = m_InputDwi->GetVectorImage()->GetLargestPossibleRegion(); - m_Parameters.m_ImageSpacing = m_InputDwi->GetVectorImage()->GetSpacing(); - m_Parameters.m_ImageOrigin = m_InputDwi->GetVectorImage()->GetOrigin(); - m_Parameters.m_ImageDirection = m_InputDwi->GetVectorImage()->GetDirection(); - m_Parameters.m_Bvalue = m_InputDwi->GetReferenceBValue(); + m_Parameters.m_SignalGen.m_ImageRegion = m_InputDwi->GetVectorImage()->GetLargestPossibleRegion(); + m_Parameters.m_SignalGen.m_ImageSpacing = m_InputDwi->GetVectorImage()->GetSpacing(); + m_Parameters.m_SignalGen.m_ImageOrigin = m_InputDwi->GetVectorImage()->GetOrigin(); + m_Parameters.m_SignalGen.m_ImageDirection = m_InputDwi->GetVectorImage()->GetDirection(); + m_Parameters.m_SignalGen.m_Bvalue = m_InputDwi->GetReferenceBValue(); m_Parameters.SetGradienDirections(m_InputDwi->GetDirections()); } bool CompareDwi(itk::VectorImage< short, 3 >* dwi1, itk::VectorImage< short, 3 >* dwi2) { typedef itk::VectorImage< short, 3 > DwiImageType; try{ itk::ImageRegionIterator< DwiImageType > it1(dwi1, dwi1->GetLargestPossibleRegion()); itk::ImageRegionIterator< DwiImageType > it2(dwi2, dwi2->GetLargestPossibleRegion()); while(!it1.IsAtEnd()) { if (it1.Get()!=it2.Get()) return false; ++it1; ++it2; } } catch(...) { return false; } return true; } void StartSimulation(string testFileName) { mitk::DiffusionImage::Pointer refImage = NULL; if (!testFileName.empty()) CPPUNIT_ASSERT(refImage = dynamic_cast*>(mitk::IOUtil::LoadDataNode(testFileName)->GetData())); itk::AddArtifactsToDwiImageFilter< short >::Pointer artifactsToDwiFilter = itk::AddArtifactsToDwiImageFilter< short >::New(); artifactsToDwiFilter->SetUseConstantRandSeed(true); artifactsToDwiFilter->SetInput(m_InputDwi->GetVectorImage()); artifactsToDwiFilter->SetParameters(m_Parameters); CPPUNIT_ASSERT_NO_THROW(artifactsToDwiFilter->Update()); mitk::DiffusionImage::Pointer testImage = mitk::DiffusionImage::New(); testImage->SetVectorImage( artifactsToDwiFilter->GetOutput() ); - testImage->SetReferenceBValue(m_Parameters.m_Bvalue); - testImage->SetDirections(m_Parameters.GetGradientDirections()); + testImage->SetReferenceBValue( m_Parameters.m_SignalGen.m_Bvalue); + testImage->SetDirections( m_Parameters.GetGradientDirections()); testImage->InitializeFromVectorImage(); if (refImage.IsNotNull()) { CPPUNIT_ASSERT_MESSAGE(testFileName, CompareDwi(testImage->GetVectorImage(), refImage->GetVectorImage())); } } void Spikes() { - m_Parameters.m_Spikes = 5; - m_Parameters.m_SpikeAmplitude = 1; + m_Parameters.m_SignalGen.m_Spikes = 5; + m_Parameters.m_SignalGen.m_SpikeAmplitude = 1; StartSimulation( GetTestDataFilePath("DiffusionImaging/Fiberfox/spikes2.dwi") ); } void GibbsRinging() { - m_Parameters.m_DoAddGibbsRinging = true; + m_Parameters.m_SignalGen.m_DoAddGibbsRinging = true; StartSimulation( GetTestDataFilePath("DiffusionImaging/Fiberfox/gibbsringing2.dwi") ); } void Ghost() { - m_Parameters.m_KspaceLineOffset = 0.25; + m_Parameters.m_SignalGen.m_KspaceLineOffset = 0.25; StartSimulation( GetTestDataFilePath("DiffusionImaging/Fiberfox/ghost2.dwi") ); } void Aliasing() { - m_Parameters.m_CroppingFactor = 0.4; + m_Parameters.m_SignalGen.m_CroppingFactor = 0.4; StartSimulation( GetTestDataFilePath("DiffusionImaging/Fiberfox/aliasing2.dwi") ); } void Eddy() { - m_Parameters.m_EddyStrength = 0.05; + m_Parameters.m_SignalGen.m_EddyStrength = 0.05; StartSimulation( GetTestDataFilePath("DiffusionImaging/Fiberfox/eddy2.dwi") ); } void RicianNoise() { mitk::RicianNoiseModel* ricianNoiseModel = new mitk::RicianNoiseModel(); ricianNoiseModel->SetNoiseVariance(1000000); ricianNoiseModel->SetSeed(0); m_Parameters.m_NoiseModel = ricianNoiseModel; StartSimulation( GetTestDataFilePath("DiffusionImaging/Fiberfox/riciannoise2.dwi") ); - delete m_Parameters.m_NoiseModel; + delete m_Parameters.m_NoiseModel; } void ChiSquareNoise() { mitk::ChiSquareNoiseModel* chiSquareNoiseModel = new mitk::ChiSquareNoiseModel(); chiSquareNoiseModel->SetNoiseVariance(1000000); chiSquareNoiseModel->SetSeed(0); m_Parameters.m_NoiseModel = chiSquareNoiseModel; StartSimulation( GetTestDataFilePath("DiffusionImaging/Fiberfox/chisquarenoise2.dwi") ); - delete m_Parameters.m_NoiseModel; + delete m_Parameters.m_NoiseModel; } void Distortions() { mitk::Image::Pointer mitkFMap = dynamic_cast(mitk::IOUtil::LoadDataNode( GetTestDataFilePath("DiffusionImaging/Fiberfox/Fieldmap.nrrd") )->GetData()); typedef itk::Image ItkDoubleImgType; ItkDoubleImgType::Pointer fMap = ItkDoubleImgType::New(); mitk::CastToItkImage(mitkFMap, fMap); - m_Parameters.m_FrequencyMap = fMap; + m_Parameters.m_SignalGen.m_FrequencyMap = fMap; StartSimulation( GetTestDataFilePath("DiffusionImaging/Fiberfox/distortions2.dwi") ); } }; MITK_TEST_SUITE_REGISTRATION(mitkFiberfoxAddArtifactsToDwi) diff --git a/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxSignalGenerationTest.cpp b/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxSignalGenerationTest.cpp index f67c4bc20c..a21db5e488 100644 --- a/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxSignalGenerationTest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxSignalGenerationTest.cpp @@ -1,276 +1,276 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include /**Documentation * Test the Fiberfox simulation functions (fiberBundle -> diffusion weighted image) */ bool CompareDwi(itk::VectorImage< short, 3 >* dwi1, itk::VectorImage< short, 3 >* dwi2) { typedef itk::VectorImage< short, 3 > DwiImageType; try{ itk::ImageRegionIterator< DwiImageType > it1(dwi1, dwi1->GetLargestPossibleRegion()); itk::ImageRegionIterator< DwiImageType > it2(dwi2, dwi2->GetLargestPossibleRegion()); while(!it1.IsAtEnd()) { if (it1.Get()!=it2.Get()) return false; ++it1; ++it2; } } catch(...) { return false; } return true; } void StartSimulation(FiberfoxParameters parameters, FiberBundleX::Pointer fiberBundle, mitk::DiffusionImage::Pointer refImage, string message) { itk::TractsToDWIImageFilter< short >::Pointer tractsToDwiFilter = itk::TractsToDWIImageFilter< short >::New(); tractsToDwiFilter->SetUseConstantRandSeed(true); tractsToDwiFilter->SetParameters(parameters); tractsToDwiFilter->SetFiberBundle(fiberBundle); tractsToDwiFilter->Update(); mitk::DiffusionImage::Pointer testImage = mitk::DiffusionImage::New(); testImage->SetVectorImage( tractsToDwiFilter->GetOutput() ); - testImage->SetReferenceBValue(parameters.m_Bvalue); + testImage->SetReferenceBValue(parameters.m_SignalGen.m_Bvalue); testImage->SetDirections(parameters.GetGradientDirections()); testImage->InitializeFromVectorImage(); if (refImage.IsNotNull()) { bool cond = CompareDwi(testImage->GetVectorImage(), refImage->GetVectorImage()); if (!cond) { MITK_INFO << "Saving test and rference image to " << mitk::IOUtil::GetTempPath(); mitk::IOUtil::SaveBaseData(testImage, mitk::IOUtil::GetTempPath()+"testImage.dwi"); mitk::IOUtil::SaveBaseData(refImage, mitk::IOUtil::GetTempPath()+"refImage.dwi"); } MITK_TEST_CONDITION_REQUIRED(cond, message); } } int mitkFiberfoxSignalGenerationTest(int argc, char* argv[]) { MITK_TEST_BEGIN("mitkFiberfoxSignalGenerationTest"); MITK_TEST_CONDITION_REQUIRED(argc>=19,"check for input data"); // input fiber bundle FiberBundleXReader::Pointer fibReader = FiberBundleXReader::New(); fibReader->SetFileName(argv[1]); fibReader->Update(); FiberBundleX::Pointer fiberBundle = dynamic_cast(fibReader->GetOutput()); // reference diffusion weighted images mitk::DiffusionImage::Pointer stickBall = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[2])->GetData()); mitk::DiffusionImage::Pointer stickAstrosticks = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[3])->GetData()); mitk::DiffusionImage::Pointer stickDot = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[4])->GetData()); mitk::DiffusionImage::Pointer tensorBall = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[5])->GetData()); mitk::DiffusionImage::Pointer stickTensorBall = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[6])->GetData()); mitk::DiffusionImage::Pointer stickTensorBallAstrosticks = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[7])->GetData()); mitk::DiffusionImage::Pointer gibbsringing = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[8])->GetData()); mitk::DiffusionImage::Pointer ghost = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[9])->GetData()); mitk::DiffusionImage::Pointer aliasing = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[10])->GetData()); mitk::DiffusionImage::Pointer eddy = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[11])->GetData()); mitk::DiffusionImage::Pointer linearmotion = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[12])->GetData()); mitk::DiffusionImage::Pointer randommotion = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[13])->GetData()); mitk::DiffusionImage::Pointer spikes = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[14])->GetData()); mitk::DiffusionImage::Pointer riciannoise = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[15])->GetData()); mitk::DiffusionImage::Pointer chisquarenoise = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[16])->GetData()); mitk::DiffusionImage::Pointer distortions = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[17])->GetData()); mitk::Image::Pointer mitkFMap = dynamic_cast(mitk::IOUtil::LoadDataNode(argv[18])->GetData()); typedef itk::Image ItkDoubleImgType; ItkDoubleImgType::Pointer fMap = ItkDoubleImgType::New(); mitk::CastToItkImage(mitkFMap, fMap); FiberfoxParameters parameters; - parameters.m_SimulateKspaceAcquisition = true; - parameters.m_SignalScale = 10000; - parameters.m_ImageRegion = stickBall->GetVectorImage()->GetLargestPossibleRegion(); - parameters.m_ImageSpacing = stickBall->GetVectorImage()->GetSpacing(); - parameters.m_ImageOrigin = stickBall->GetVectorImage()->GetOrigin(); - parameters.m_ImageDirection = stickBall->GetVectorImage()->GetDirection(); - parameters.m_Bvalue = stickBall->GetReferenceBValue(); + parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; + parameters.m_SignalGen.m_SignalScale = 10000; + parameters.m_SignalGen.m_ImageRegion = stickBall->GetVectorImage()->GetLargestPossibleRegion(); + parameters.m_SignalGen.m_ImageSpacing = stickBall->GetVectorImage()->GetSpacing(); + parameters.m_SignalGen.m_ImageOrigin = stickBall->GetVectorImage()->GetOrigin(); + parameters.m_SignalGen.m_ImageDirection = stickBall->GetVectorImage()->GetDirection(); + parameters.m_SignalGen.m_Bvalue = stickBall->GetReferenceBValue(); parameters.SetGradienDirections(stickBall->GetDirections()); // intra and inter axonal compartments mitk::StickModel stickModel; - stickModel.SetBvalue(parameters.m_Bvalue); + stickModel.SetBvalue(parameters.m_SignalGen.m_Bvalue); stickModel.SetT2(110); stickModel.SetDiffusivity(0.001); stickModel.SetGradientList(parameters.GetGradientDirections()); mitk::TensorModel tensorModel; tensorModel.SetT2(110); - stickModel.SetBvalue(parameters.m_Bvalue); + stickModel.SetBvalue(parameters.m_SignalGen.m_Bvalue); tensorModel.SetDiffusivity1(0.001); tensorModel.SetDiffusivity2(0.00025); tensorModel.SetDiffusivity3(0.00025); tensorModel.SetGradientList(parameters.GetGradientDirections()); // extra axonal compartment models mitk::BallModel ballModel; ballModel.SetT2(80); - ballModel.SetBvalue(parameters.m_Bvalue); + ballModel.SetBvalue(parameters.m_SignalGen.m_Bvalue); ballModel.SetDiffusivity(0.001); ballModel.SetGradientList(parameters.GetGradientDirections()); mitk::AstroStickModel astrosticksModel; astrosticksModel.SetT2(80); - astrosticksModel.SetBvalue(parameters.m_Bvalue); + astrosticksModel.SetBvalue(parameters.m_SignalGen.m_Bvalue); astrosticksModel.SetDiffusivity(0.001); astrosticksModel.SetRandomizeSticks(true); astrosticksModel.SetSeed(0); astrosticksModel.SetGradientList(parameters.GetGradientDirections()); mitk::DotModel dotModel; dotModel.SetT2(80); dotModel.SetGradientList(parameters.GetGradientDirections()); // noise models mitk::RicianNoiseModel* ricianNoiseModel = new mitk::RicianNoiseModel(); ricianNoiseModel->SetNoiseVariance(1000000); ricianNoiseModel->SetSeed(0); // Rician noise mitk::ChiSquareNoiseModel* chiSquareNoiseModel = new mitk::ChiSquareNoiseModel(); chiSquareNoiseModel->SetNoiseVariance(1000000); chiSquareNoiseModel->SetSeed(0); try{ // Stick-Ball parameters.m_FiberModelList.push_back(&stickModel); parameters.m_NonFiberModelList.push_back(&ballModel); StartSimulation(parameters, fiberBundle, stickBall, argv[2]); // Srick-Astrosticks parameters.m_NonFiberModelList.clear(); parameters.m_NonFiberModelList.push_back(&astrosticksModel); StartSimulation(parameters, fiberBundle, stickAstrosticks, argv[3]); // Stick-Dot parameters.m_NonFiberModelList.clear(); parameters.m_NonFiberModelList.push_back(&dotModel); StartSimulation(parameters, fiberBundle, stickDot, argv[4]); // Tensor-Ball parameters.m_FiberModelList.clear(); parameters.m_FiberModelList.push_back(&tensorModel); parameters.m_NonFiberModelList.clear(); parameters.m_NonFiberModelList.push_back(&ballModel); StartSimulation(parameters, fiberBundle, tensorBall, argv[5]); // Stick-Tensor-Ball parameters.m_FiberModelList.clear(); parameters.m_FiberModelList.push_back(&stickModel); parameters.m_FiberModelList.push_back(&tensorModel); parameters.m_NonFiberModelList.clear(); parameters.m_NonFiberModelList.push_back(&ballModel); StartSimulation(parameters, fiberBundle, stickTensorBall, argv[6]); // Stick-Tensor-Ball-Astrosticks parameters.m_NonFiberModelList.push_back(&astrosticksModel); StartSimulation(parameters, fiberBundle, stickTensorBallAstrosticks, argv[7]); // Gibbs ringing parameters.m_FiberModelList.clear(); parameters.m_FiberModelList.push_back(&stickModel); parameters.m_NonFiberModelList.clear(); parameters.m_NonFiberModelList.push_back(&ballModel); - parameters.m_DoAddGibbsRinging = true; + parameters.m_SignalGen.m_DoAddGibbsRinging = true; StartSimulation(parameters, fiberBundle, gibbsringing, argv[8]); // Ghost - parameters.m_DoAddGibbsRinging = false; - parameters.m_KspaceLineOffset = 0.25; + parameters.m_SignalGen.m_DoAddGibbsRinging = false; + parameters.m_SignalGen.m_KspaceLineOffset = 0.25; StartSimulation(parameters, fiberBundle, ghost, argv[9]); // Aliasing - parameters.m_KspaceLineOffset = 0; - parameters.m_CroppingFactor = 0.4; - parameters.m_SignalScale = 1000; + parameters.m_SignalGen.m_KspaceLineOffset = 0; + parameters.m_SignalGen.m_CroppingFactor = 0.4; + parameters.m_SignalGen.m_SignalScale = 1000; StartSimulation(parameters, fiberBundle, aliasing, argv[10]); // Eddy currents - parameters.m_CroppingFactor = 1; - parameters.m_SignalScale = 10000; - parameters.m_EddyStrength = 0.05; + parameters.m_SignalGen.m_CroppingFactor = 1; + parameters.m_SignalGen.m_SignalScale = 10000; + parameters.m_SignalGen.m_EddyStrength = 0.05; StartSimulation(parameters, fiberBundle, eddy, argv[11]); // Motion (linear) - parameters.m_EddyStrength = 0.0; - parameters.m_DoAddMotion = true; - parameters.m_DoRandomizeMotion = false; - parameters.m_Translation[1] = 10; - parameters.m_Rotation[2] = 90; + parameters.m_SignalGen.m_EddyStrength = 0.0; + parameters.m_SignalGen.m_DoAddMotion = true; + parameters.m_SignalGen.m_DoRandomizeMotion = false; + parameters.m_SignalGen.m_Translation[1] = 10; + parameters.m_SignalGen.m_Rotation[2] = 90; StartSimulation(parameters, fiberBundle, linearmotion, argv[12]); // Motion (random) - parameters.m_DoRandomizeMotion = true; - parameters.m_Translation[1] = 5; - parameters.m_Rotation[2] = 45; + parameters.m_SignalGen.m_DoRandomizeMotion = true; + parameters.m_SignalGen.m_Translation[1] = 5; + parameters.m_SignalGen.m_Rotation[2] = 45; StartSimulation(parameters, fiberBundle, randommotion, argv[13]); // Spikes - parameters.m_DoAddMotion = false; - parameters.m_Spikes = 5; - parameters.m_SpikeAmplitude = 1; + parameters.m_SignalGen.m_DoAddMotion = false; + parameters.m_SignalGen.m_Spikes = 5; + parameters.m_SignalGen.m_SpikeAmplitude = 1; StartSimulation(parameters, fiberBundle, spikes, argv[14]); // Rician noise - parameters.m_Spikes = 0; + parameters.m_SignalGen.m_Spikes = 0; parameters.m_NoiseModel = ricianNoiseModel; StartSimulation(parameters, fiberBundle, riciannoise, argv[15]); delete parameters.m_NoiseModel; // Chi-square noise parameters.m_NoiseModel = chiSquareNoiseModel; StartSimulation(parameters, fiberBundle, chisquarenoise, argv[16]); delete parameters.m_NoiseModel; // Distortions parameters.m_NoiseModel = NULL; - parameters.m_FrequencyMap = fMap; + parameters.m_SignalGen.m_FrequencyMap = fMap; StartSimulation(parameters, fiberBundle, distortions, argv[17]); } catch (std::exception &e) { MITK_TEST_CONDITION_REQUIRED(false, e.what()); } // always end with this! MITK_TEST_END(); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.cpp index cbc63205ee..222ede5bf9 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.cpp @@ -1,2603 +1,2604 @@ /*=================================================================== 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. ===================================================================*/ //misc #define _USE_MATH_DEFINES #include // Blueberry #include #include // Qmitk #include "QmitkFiberfoxView.h" // MITK #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "usModuleRegistry.h" #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include QmitkFiberfoxWorker::QmitkFiberfoxWorker(QmitkFiberfoxView* view) : m_View(view) { } void QmitkFiberfoxWorker::run() { try{ switch (m_FilterType) { case 0: m_View->m_TractsToDwiFilter->Update(); break; case 1: m_View->m_ArtifactsToDwiFilter->Update(); break; } } catch( ... ) { } m_View->m_Thread.quit(); } const std::string QmitkFiberfoxView::VIEW_ID = "org.mitk.views.fiberfoxview"; QmitkFiberfoxView::QmitkFiberfoxView() : QmitkAbstractView() , m_Controls( 0 ) , m_SelectedImage( NULL ) , m_Worker(this) , m_ThreadIsRunning(false) { m_Worker.moveToThread(&m_Thread); connect(&m_Thread, SIGNAL(started()), this, SLOT(BeforeThread())); connect(&m_Thread, SIGNAL(started()), &m_Worker, SLOT(run())); connect(&m_Thread, SIGNAL(finished()), this, SLOT(AfterThread())); connect(&m_Thread, SIGNAL(terminated()), this, SLOT(AfterThread())); m_SimulationTimer = new QTimer(this); m_StickModel1.m_CompartmentId = 1; m_StickModel2.m_CompartmentId = 2; m_TensorModel1.m_CompartmentId = 1; m_TensorModel2.m_CompartmentId = 2; m_ZeppelinModel1.m_CompartmentId = 1; m_ZeppelinModel2.m_CompartmentId = 2; m_BallModel1.m_CompartmentId = 3; m_BallModel2.m_CompartmentId = 4; m_AstrosticksModel1.m_CompartmentId = 3; m_AstrosticksModel2.m_CompartmentId = 4; m_DotModel1.m_CompartmentId = 3; m_DotModel2.m_CompartmentId = 4; m_PrototypeModel1.m_CompartmentId = 1; m_PrototypeModel3.m_CompartmentId = 3; m_PrototypeModel4.m_CompartmentId = 4; } void QmitkFiberfoxView::KillThread() { MITK_INFO << "Aborting DWI simulation."; switch (m_Worker.m_FilterType) { case 0: m_TractsToDwiFilter->SetAbortGenerateData(true); break; case 1: m_ArtifactsToDwiFilter->SetAbortGenerateData(true); break; } m_Controls->m_AbortSimulationButton->setEnabled(false); m_Controls->m_AbortSimulationButton->setText("Aborting simulation ..."); } void QmitkFiberfoxView::BeforeThread() { m_SimulationTime = QTime::currentTime(); m_SimulationTimer->start(100); m_Controls->m_AbortSimulationButton->setVisible(true); m_Controls->m_GenerateImageButton->setVisible(false); m_Controls->m_SimulationStatusText->setVisible(true); m_ThreadIsRunning = true; } void QmitkFiberfoxView::AfterThread() { UpdateSimulationStatus(); m_SimulationTimer->stop(); m_Controls->m_AbortSimulationButton->setVisible(false); m_Controls->m_AbortSimulationButton->setEnabled(true); m_Controls->m_AbortSimulationButton->setText("Abort simulation"); m_Controls->m_GenerateImageButton->setVisible(true); m_ThreadIsRunning = false; QString statusText; FiberfoxParameters parameters; mitk::DiffusionImage::Pointer mitkImage = mitk::DiffusionImage::New(); switch (m_Worker.m_FilterType) { case 0: { statusText = QString(m_TractsToDwiFilter->GetStatusText().c_str()); if (m_TractsToDwiFilter->GetAbortGenerateData()) { MITK_INFO << "Simulation aborted."; return; } parameters = m_TractsToDwiFilter->GetParameters(); mitkImage->SetVectorImage( m_TractsToDwiFilter->GetOutput() ); - mitkImage->SetReferenceBValue(parameters.m_Bvalue); + mitkImage->SetReferenceBValue(parameters.m_SignalGen.m_Bvalue); mitkImage->SetDirections(parameters.GetGradientDirections()); mitkImage->InitializeFromVectorImage(); - parameters.m_ResultNode->SetData( mitkImage ); + parameters.m_Misc.m_ResultNode->SetData( mitkImage ); - parameters.m_ResultNode->SetName(parameters.m_ParentNode->GetName() - +"_D"+QString::number(parameters.m_ImageRegion.GetSize(0)).toStdString() - +"-"+QString::number(parameters.m_ImageRegion.GetSize(1)).toStdString() - +"-"+QString::number(parameters.m_ImageRegion.GetSize(2)).toStdString() - +"_S"+QString::number(parameters.m_ImageSpacing[0]).toStdString() - +"-"+QString::number(parameters.m_ImageSpacing[1]).toStdString() - +"-"+QString::number(parameters.m_ImageSpacing[2]).toStdString() - +"_b"+QString::number(parameters.m_Bvalue).toStdString() - +"_"+parameters.m_SignalModelString - +parameters.m_ArtifactModelString); + parameters.m_Misc.m_ResultNode->SetName(parameters.m_Misc.m_ParentNode->GetName() + +"_D"+QString::number(parameters.m_SignalGen.m_ImageRegion.GetSize(0)).toStdString() + +"-"+QString::number(parameters.m_SignalGen.m_ImageRegion.GetSize(1)).toStdString() + +"-"+QString::number(parameters.m_SignalGen.m_ImageRegion.GetSize(2)).toStdString() + +"_S"+QString::number(parameters.m_SignalGen.m_ImageSpacing[0]).toStdString() + +"-"+QString::number(parameters.m_SignalGen.m_ImageSpacing[1]).toStdString() + +"-"+QString::number(parameters.m_SignalGen.m_ImageSpacing[2]).toStdString() + +"_b"+QString::number(parameters.m_SignalGen.m_Bvalue).toStdString() + +"_"+parameters.m_Misc.m_SignalModelString + +parameters.m_Misc.m_ArtifactModelString); - GetDataStorage()->Add(parameters.m_ResultNode, parameters.m_ParentNode); + GetDataStorage()->Add(parameters.m_Misc.m_ResultNode, parameters.m_Misc.m_ParentNode); - parameters.m_ResultNode->SetProperty( "levelwindow", mitk::LevelWindowProperty::New(m_TractsToDwiFilter->GetLevelWindow()) ); + parameters.m_Misc.m_ResultNode->SetProperty( "levelwindow", mitk::LevelWindowProperty::New(m_TractsToDwiFilter->GetLevelWindow()) ); if (m_Controls->m_VolumeFractionsBox->isChecked()) { std::vector< itk::TractsToDWIImageFilter< short >::ItkDoubleImgType::Pointer > volumeFractions = m_TractsToDwiFilter->GetVolumeFractions(); for (unsigned int k=0; kInitializeByItk(volumeFractions.at(k).GetPointer()); image->SetVolume(volumeFractions.at(k)->GetBufferPointer()); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); - node->SetName(parameters.m_ParentNode->GetName()+"_CompartmentVolume-"+QString::number(k).toStdString()); - GetDataStorage()->Add(node, parameters.m_ParentNode); + node->SetName(parameters.m_Misc.m_ParentNode->GetName()+"_CompartmentVolume-"+QString::number(k).toStdString()); + GetDataStorage()->Add(node, parameters.m_Misc.m_ParentNode); } } m_TractsToDwiFilter = NULL; break; } case 1: { statusText = QString(m_ArtifactsToDwiFilter->GetStatusText().c_str()); if (m_ArtifactsToDwiFilter->GetAbortGenerateData()) { MITK_INFO << "Simulation aborted."; return; } parameters = m_ArtifactsToDwiFilter->GetParameters().CopyParameters(); - mitk::DiffusionImage::Pointer diffImg = dynamic_cast*>(parameters.m_ParentNode->GetData()); + mitk::DiffusionImage::Pointer diffImg = dynamic_cast*>(parameters.m_Misc.m_ParentNode->GetData()); mitkImage = mitk::DiffusionImage::New(); mitkImage->SetVectorImage( m_ArtifactsToDwiFilter->GetOutput() ); mitkImage->SetReferenceBValue(diffImg->GetReferenceBValue()); mitkImage->SetDirections(diffImg->GetDirections()); mitkImage->InitializeFromVectorImage(); - parameters.m_ResultNode->SetData( mitkImage ); - parameters.m_ResultNode->SetName(parameters.m_ParentNode->GetName()+parameters.m_ArtifactModelString); - GetDataStorage()->Add(parameters.m_ResultNode, parameters.m_ParentNode); + parameters.m_Misc.m_ResultNode->SetData( mitkImage ); + parameters.m_Misc.m_ResultNode->SetName(parameters.m_Misc.m_ParentNode->GetName()+parameters.m_Misc.m_ArtifactModelString); + GetDataStorage()->Add(parameters.m_Misc.m_ResultNode, parameters.m_Misc.m_ParentNode); m_ArtifactsToDwiFilter = NULL; break; } } - mitk::BaseData::Pointer basedata = parameters.m_ResultNode->GetData(); + 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_OutputPath.empty()) + if (!parameters.m_Misc.m_OutputPath.empty()) { try{ - QString outputFileName(parameters.m_OutputPath.c_str()); - outputFileName += parameters.m_ResultNode->GetName().c_str(); + QString outputFileName(parameters.m_Misc.m_OutputPath.c_str()); + outputFileName += parameters.m_Misc.m_ResultNode->GetName().c_str(); outputFileName.replace(QString("."), QString("_")); outputFileName += ".dwi"; QString status("Saving output image to "); status += outputFileName; m_Controls->m_SimulationStatusText->append(status); mitk::IOUtil::SaveBaseData(mitkImage, outputFileName.toStdString()); m_Controls->m_SimulationStatusText->append("File saved successfully."); } catch (itk::ExceptionObject &e) { QString status("Exception during DWI writing: "); status += e.GetDescription(); m_Controls->m_SimulationStatusText->append(status); } catch (...) { m_Controls->m_SimulationStatusText->append("Unknown exception during DWI writing!"); } } - parameters.m_FrequencyMap = NULL; + parameters.m_SignalGen.m_FrequencyMap = NULL; } void QmitkFiberfoxView::UpdateSimulationStatus() { QString statusText; switch (m_Worker.m_FilterType) { case 0: statusText = QString(m_TractsToDwiFilter->GetStatusText().c_str()); break; case 1: statusText = QString(m_ArtifactsToDwiFilter->GetStatusText().c_str()); break; } if (QString::compare(m_SimulationStatusText,statusText)!=0) { m_Controls->m_SimulationStatusText->clear(); statusText = "
"+statusText+"
"; m_Controls->m_SimulationStatusText->setText(statusText); QScrollBar *vScrollBar = m_Controls->m_SimulationStatusText->verticalScrollBar(); vScrollBar->triggerAction(QScrollBar::SliderToMaximum); } } // Destructor QmitkFiberfoxView::~QmitkFiberfoxView() { delete m_SimulationTimer; } void QmitkFiberfoxView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkFiberfoxViewControls; m_Controls->setupUi( parent ); m_Controls->m_StickWidget1->setVisible(true); m_Controls->m_StickWidget2->setVisible(false); m_Controls->m_ZeppelinWidget1->setVisible(false); m_Controls->m_ZeppelinWidget2->setVisible(false); m_Controls->m_TensorWidget1->setVisible(false); m_Controls->m_TensorWidget2->setVisible(false); m_Controls->m_BallWidget1->setVisible(true); m_Controls->m_BallWidget2->setVisible(false); m_Controls->m_AstrosticksWidget1->setVisible(false); m_Controls->m_AstrosticksWidget2->setVisible(false); m_Controls->m_DotWidget1->setVisible(false); m_Controls->m_DotWidget2->setVisible(false); m_Controls->m_PrototypeWidget1->setVisible(false); m_Controls->m_PrototypeWidget2->setVisible(false); m_Controls->m_PrototypeWidget3->setVisible(false); m_Controls->m_PrototypeWidget4->setVisible(false); m_Controls->m_PrototypeWidget3->SetMinFa(0.0); m_Controls->m_PrototypeWidget3->SetMaxFa(0.15); m_Controls->m_PrototypeWidget4->SetMinFa(0.0); m_Controls->m_PrototypeWidget4->SetMaxFa(0.15); m_Controls->m_PrototypeWidget3->SetMinAdc(0.0); m_Controls->m_PrototypeWidget3->SetMaxAdc(0.001); m_Controls->m_PrototypeWidget4->SetMinAdc(0.003); m_Controls->m_PrototypeWidget4->SetMaxAdc(0.004); m_Controls->m_Comp4FractionFrame->setVisible(false); m_Controls->m_DiffusionPropsMessage->setVisible(false); m_Controls->m_GeometryMessage->setVisible(false); m_Controls->m_AdvancedSignalOptionsFrame->setVisible(false); m_Controls->m_AdvancedFiberOptionsFrame->setVisible(false); m_Controls->m_VarianceBox->setVisible(false); m_Controls->m_NoiseFrame->setVisible(false); m_Controls->m_GhostFrame->setVisible(false); m_Controls->m_DistortionsFrame->setVisible(false); m_Controls->m_EddyFrame->setVisible(false); m_Controls->m_SpikeFrame->setVisible(false); m_Controls->m_AliasingFrame->setVisible(false); m_Controls->m_MotionArtifactFrame->setVisible(false); m_ParameterFile = QDir::currentPath()+"/param.ffp"; m_Controls->m_AbortSimulationButton->setVisible(false); m_Controls->m_SimulationStatusText->setVisible(false); m_Controls->m_FrequencyMapBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isMitkImage = mitk::TNodePredicateDataType::New(); mitk::NodePredicateDataType::Pointer isDwi = mitk::NodePredicateDataType::New("DiffusionImage"); mitk::NodePredicateDataType::Pointer isDti = mitk::NodePredicateDataType::New("TensorImage"); mitk::NodePredicateDataType::Pointer isQbi = mitk::NodePredicateDataType::New("QBallImage"); mitk::NodePredicateOr::Pointer isDiffusionImage = mitk::NodePredicateOr::New(isDwi, isDti); isDiffusionImage = mitk::NodePredicateOr::New(isDiffusionImage, isQbi); mitk::NodePredicateNot::Pointer noDiffusionImage = mitk::NodePredicateNot::New(isDiffusionImage); mitk::NodePredicateAnd::Pointer finalPredicate = mitk::NodePredicateAnd::New(isMitkImage, noDiffusionImage); m_Controls->m_FrequencyMapBox->SetPredicate(finalPredicate); m_Controls->m_Comp4VolumeFraction->SetDataStorage(this->GetDataStorage()); m_Controls->m_Comp4VolumeFraction->SetPredicate(finalPredicate); connect( m_SimulationTimer, SIGNAL(timeout()), this, SLOT(UpdateSimulationStatus()) ); connect((QObject*) m_Controls->m_AbortSimulationButton, SIGNAL(clicked()), (QObject*) this, SLOT(KillThread())); connect((QObject*) m_Controls->m_GenerateImageButton, SIGNAL(clicked()), (QObject*) this, SLOT(GenerateImage())); connect((QObject*) m_Controls->m_GenerateFibersButton, SIGNAL(clicked()), (QObject*) this, SLOT(GenerateFibers())); connect((QObject*) m_Controls->m_CircleButton, SIGNAL(clicked()), (QObject*) this, SLOT(OnDrawROI())); connect((QObject*) m_Controls->m_FlipButton, SIGNAL(clicked()), (QObject*) this, SLOT(OnFlipButton())); connect((QObject*) m_Controls->m_JoinBundlesButton, SIGNAL(clicked()), (QObject*) this, SLOT(JoinBundles())); connect((QObject*) m_Controls->m_VarianceBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnVarianceChanged(double))); connect((QObject*) m_Controls->m_DistributionBox, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(OnDistributionChanged(int))); connect((QObject*) m_Controls->m_FiberDensityBox, SIGNAL(valueChanged(int)), (QObject*) this, SLOT(OnFiberDensityChanged(int))); connect((QObject*) m_Controls->m_FiberSamplingBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnFiberSamplingChanged(double))); connect((QObject*) m_Controls->m_TensionBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnTensionChanged(double))); connect((QObject*) m_Controls->m_ContinuityBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnContinuityChanged(double))); connect((QObject*) m_Controls->m_BiasBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnBiasChanged(double))); connect((QObject*) m_Controls->m_AddNoise, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddNoise(int))); connect((QObject*) m_Controls->m_AddGhosts, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddGhosts(int))); connect((QObject*) m_Controls->m_AddDistortions, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddDistortions(int))); connect((QObject*) m_Controls->m_AddEddy, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddEddy(int))); connect((QObject*) m_Controls->m_AddSpikes, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddSpikes(int))); connect((QObject*) m_Controls->m_AddAliasing, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddAliasing(int))); connect((QObject*) m_Controls->m_AddMotion, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddMotion(int))); connect((QObject*) m_Controls->m_ConstantRadiusBox, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnConstantRadius(int))); connect((QObject*) m_Controls->m_CopyBundlesButton, SIGNAL(clicked()), (QObject*) this, SLOT(CopyBundles())); connect((QObject*) m_Controls->m_TransformBundlesButton, SIGNAL(clicked()), (QObject*) this, SLOT(ApplyTransform())); connect((QObject*) m_Controls->m_AlignOnGrid, SIGNAL(clicked()), (QObject*) this, SLOT(AlignOnGrid())); connect((QObject*) m_Controls->m_Compartment1Box, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(Comp1ModelFrameVisibility(int))); connect((QObject*) m_Controls->m_Compartment2Box, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(Comp2ModelFrameVisibility(int))); connect((QObject*) m_Controls->m_Compartment3Box, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(Comp3ModelFrameVisibility(int))); connect((QObject*) m_Controls->m_Compartment4Box, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(Comp4ModelFrameVisibility(int))); connect((QObject*) m_Controls->m_AdvancedOptionsBox, SIGNAL( stateChanged(int)), (QObject*) this, SLOT(ShowAdvancedOptions(int))); connect((QObject*) m_Controls->m_AdvancedOptionsBox_2, SIGNAL( stateChanged(int)), (QObject*) this, SLOT(ShowAdvancedOptions(int))); connect((QObject*) m_Controls->m_SaveParametersButton, SIGNAL(clicked()), (QObject*) this, SLOT(SaveParameters())); connect((QObject*) m_Controls->m_LoadParametersButton, SIGNAL(clicked()), (QObject*) this, SLOT(LoadParameters())); connect((QObject*) m_Controls->m_OutputPathButton, SIGNAL(clicked()), (QObject*) this, SLOT(SetOutputPath())); } } template< class ScalarType > FiberfoxParameters< ScalarType > QmitkFiberfoxView::UpdateImageParameters() { FiberfoxParameters< ScalarType > parameters; - parameters.m_OutputPath = ""; - parameters.m_FiberGenerationParameters.m_AdvancedOptions = m_Controls->m_AdvancedOptionsBox->isChecked(); - parameters.m_AdvancedOptions = m_Controls->m_AdvancedOptionsBox_2->isChecked(); + parameters.m_Misc.m_OutputPath = ""; + parameters.m_FiberGen.m_AdvancedOptions = m_Controls->m_AdvancedOptionsBox->isChecked(); + parameters.m_Misc.m_AdvancedOptions = m_Controls->m_AdvancedOptionsBox_2->isChecked(); string outputPath = m_Controls->m_SavePathEdit->text().toStdString(); if (outputPath.compare("-")!=0) { - parameters.m_OutputPath = outputPath; - parameters.m_OutputPath += "/"; + parameters.m_Misc.m_OutputPath = outputPath; + parameters.m_Misc.m_OutputPath += "/"; } if (m_MaskImageNode.IsNotNull()) { mitk::Image::Pointer mitkMaskImage = dynamic_cast(m_MaskImageNode->GetData()); - mitk::CastToItkImage(mitkMaskImage, parameters.m_MaskImage); + mitk::CastToItkImage(mitkMaskImage, parameters.m_SignalGen.m_MaskImage); itk::ImageDuplicator::Pointer duplicator = itk::ImageDuplicator::New(); - duplicator->SetInputImage(parameters.m_MaskImage); + duplicator->SetInputImage(parameters.m_SignalGen.m_MaskImage); duplicator->Update(); - parameters.m_MaskImage = duplicator->GetOutput(); + parameters.m_SignalGen.m_MaskImage = duplicator->GetOutput(); } if (m_SelectedDWI.IsNotNull()) // use parameters of selected DWI { mitk::DiffusionImage::Pointer dwi = dynamic_cast*>(m_SelectedDWI->GetData()); - parameters.m_ImageRegion = dwi->GetVectorImage()->GetLargestPossibleRegion(); - parameters.m_ImageSpacing = dwi->GetVectorImage()->GetSpacing(); - parameters.m_ImageOrigin = dwi->GetVectorImage()->GetOrigin(); - parameters.m_ImageDirection = dwi->GetVectorImage()->GetDirection(); - parameters.m_Bvalue = dwi->GetReferenceBValue(); + parameters.m_SignalGen.m_ImageRegion = dwi->GetVectorImage()->GetLargestPossibleRegion(); + parameters.m_SignalGen.m_ImageSpacing = dwi->GetVectorImage()->GetSpacing(); + parameters.m_SignalGen.m_ImageOrigin = dwi->GetVectorImage()->GetOrigin(); + parameters.m_SignalGen.m_ImageDirection = dwi->GetVectorImage()->GetDirection(); + parameters.m_SignalGen.m_Bvalue = dwi->GetReferenceBValue(); parameters.SetGradienDirections(dwi->GetDirections()); } else if (m_SelectedImage.IsNotNull()) // use geometry of selected image { mitk::Image::Pointer img = dynamic_cast(m_SelectedImage->GetData()); itk::Image< float, 3 >::Pointer itkImg = itk::Image< float, 3 >::New(); CastToItkImage< itk::Image< float, 3 > >(img, itkImg); - parameters.m_ImageRegion = itkImg->GetLargestPossibleRegion(); - parameters.m_ImageSpacing = itkImg->GetSpacing(); - parameters.m_ImageOrigin = itkImg->GetOrigin(); - parameters.m_ImageDirection = itkImg->GetDirection(); + parameters.m_SignalGen.m_ImageRegion = itkImg->GetLargestPossibleRegion(); + parameters.m_SignalGen.m_ImageSpacing = itkImg->GetSpacing(); + parameters.m_SignalGen.m_ImageOrigin = itkImg->GetOrigin(); + parameters.m_SignalGen.m_ImageDirection = itkImg->GetDirection(); parameters.SetNumWeightedGradients(m_Controls->m_NumGradientsBox->value()); - parameters.m_Bvalue = m_Controls->m_BvalueBox->value(); + parameters.m_SignalGen.m_Bvalue = m_Controls->m_BvalueBox->value(); } else // use GUI parameters { - parameters.m_ImageRegion.SetSize(0, m_Controls->m_SizeX->value()); - parameters.m_ImageRegion.SetSize(1, m_Controls->m_SizeY->value()); - parameters.m_ImageRegion.SetSize(2, m_Controls->m_SizeZ->value()); - parameters.m_ImageSpacing[0] = m_Controls->m_SpacingX->value(); - parameters.m_ImageSpacing[1] = m_Controls->m_SpacingY->value(); - parameters.m_ImageSpacing[2] = m_Controls->m_SpacingZ->value(); - parameters.m_ImageOrigin[0] = parameters.m_ImageSpacing[0]/2; - parameters.m_ImageOrigin[1] = parameters.m_ImageSpacing[1]/2; - parameters.m_ImageOrigin[2] = parameters.m_ImageSpacing[2]/2; - parameters.m_ImageDirection.SetIdentity(); + parameters.m_SignalGen.m_ImageRegion.SetSize(0, m_Controls->m_SizeX->value()); + parameters.m_SignalGen.m_ImageRegion.SetSize(1, m_Controls->m_SizeY->value()); + parameters.m_SignalGen.m_ImageRegion.SetSize(2, m_Controls->m_SizeZ->value()); + parameters.m_SignalGen.m_ImageSpacing[0] = m_Controls->m_SpacingX->value(); + parameters.m_SignalGen.m_ImageSpacing[1] = m_Controls->m_SpacingY->value(); + parameters.m_SignalGen.m_ImageSpacing[2] = m_Controls->m_SpacingZ->value(); + parameters.m_SignalGen.m_ImageOrigin[0] = parameters.m_SignalGen.m_ImageSpacing[0]/2; + parameters.m_SignalGen.m_ImageOrigin[1] = parameters.m_SignalGen.m_ImageSpacing[1]/2; + parameters.m_SignalGen.m_ImageOrigin[2] = parameters.m_SignalGen.m_ImageSpacing[2]/2; + parameters.m_SignalGen.m_ImageDirection.SetIdentity(); parameters.SetNumWeightedGradients(m_Controls->m_NumGradientsBox->value()); - parameters.m_Bvalue = m_Controls->m_BvalueBox->value(); + parameters.m_SignalGen.m_Bvalue = m_Controls->m_BvalueBox->value(); parameters.GenerateGradientHalfShell(); } // signal relaxation - parameters.m_DoSimulateRelaxation = m_Controls->m_RelaxationBox->isChecked(); - parameters.m_SimulateKspaceAcquisition = parameters.m_DoSimulateRelaxation; - if (parameters.m_DoSimulateRelaxation && m_SelectedBundles.size()>0 ) - parameters.m_ArtifactModelString += "_RELAX"; + parameters.m_SignalGen.m_DoSimulateRelaxation = m_Controls->m_RelaxationBox->isChecked(); + parameters.m_SignalGen.m_SimulateKspaceAcquisition = parameters.m_SignalGen.m_DoSimulateRelaxation; + if (parameters.m_SignalGen.m_DoSimulateRelaxation && m_SelectedBundles.size()>0 ) + parameters.m_Misc.m_ArtifactModelString += "_RELAX"; // N/2 ghosts - parameters.m_DoAddGhosts = m_Controls->m_AddGhosts->isChecked(); + parameters.m_SignalGen.m_DoAddGhosts = m_Controls->m_AddGhosts->isChecked(); if (m_Controls->m_AddGhosts->isChecked()) { - parameters.m_SimulateKspaceAcquisition = true; - parameters.m_ArtifactModelString += "_GHOST"; - parameters.m_KspaceLineOffset = m_Controls->m_kOffsetBox->value(); - parameters.m_ResultNode->AddProperty("Fiberfox.Ghost", DoubleProperty::New(parameters.m_KspaceLineOffset)); + parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; + parameters.m_Misc.m_ArtifactModelString += "_GHOST"; + parameters.m_SignalGen.m_KspaceLineOffset = m_Controls->m_kOffsetBox->value(); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Ghost", DoubleProperty::New(parameters.m_SignalGen.m_KspaceLineOffset)); } else - parameters.m_KspaceLineOffset = 0; + parameters.m_SignalGen.m_KspaceLineOffset = 0; // Aliasing - parameters.m_DoAddAliasing = m_Controls->m_AddAliasing->isChecked(); + parameters.m_SignalGen.m_DoAddAliasing = m_Controls->m_AddAliasing->isChecked(); if (m_Controls->m_AddAliasing->isChecked()) { - parameters.m_SimulateKspaceAcquisition = true; - parameters.m_ArtifactModelString += "_ALIASING"; - parameters.m_CroppingFactor = (100-m_Controls->m_WrapBox->value())/100; - parameters.m_ResultNode->AddProperty("Fiberfox.Aliasing", DoubleProperty::New(m_Controls->m_WrapBox->value())); + parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; + parameters.m_Misc.m_ArtifactModelString += "_ALIASING"; + parameters.m_SignalGen.m_CroppingFactor = (100-m_Controls->m_WrapBox->value())/100; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Aliasing", DoubleProperty::New(m_Controls->m_WrapBox->value())); } // Spikes - parameters.m_DoAddSpikes = m_Controls->m_AddSpikes->isChecked(); + parameters.m_SignalGen.m_DoAddSpikes = m_Controls->m_AddSpikes->isChecked(); if (m_Controls->m_AddSpikes->isChecked()) { - parameters.m_SimulateKspaceAcquisition = true; - parameters.m_Spikes = m_Controls->m_SpikeNumBox->value(); - parameters.m_SpikeAmplitude = m_Controls->m_SpikeScaleBox->value(); - parameters.m_ArtifactModelString += "_SPIKES"; - parameters.m_ResultNode->AddProperty("Fiberfox.Spikes.Number", IntProperty::New(parameters.m_Spikes)); - parameters.m_ResultNode->AddProperty("Fiberfox.Spikes.Amplitude", DoubleProperty::New(parameters.m_SpikeAmplitude)); + parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; + parameters.m_SignalGen.m_Spikes = m_Controls->m_SpikeNumBox->value(); + parameters.m_SignalGen.m_SpikeAmplitude = m_Controls->m_SpikeScaleBox->value(); + parameters.m_Misc.m_ArtifactModelString += "_SPIKES"; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Spikes.Number", IntProperty::New(parameters.m_SignalGen.m_Spikes)); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Spikes.Amplitude", DoubleProperty::New(parameters.m_SignalGen.m_SpikeAmplitude)); } // gibbs ringing - parameters.m_DoAddGibbsRinging = m_Controls->m_AddGibbsRinging->isChecked(); + parameters.m_SignalGen.m_DoAddGibbsRinging = m_Controls->m_AddGibbsRinging->isChecked(); if (m_Controls->m_AddGibbsRinging->isChecked()) { - parameters.m_SimulateKspaceAcquisition = true; - parameters.m_ResultNode->AddProperty("Fiberfox.Ringing", BoolProperty::New(true)); - parameters.m_ArtifactModelString += "_RINGING"; + parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Ringing", BoolProperty::New(true)); + parameters.m_Misc.m_ArtifactModelString += "_RINGING"; } // add distortions - parameters.m_DoAddDistortions = m_Controls->m_AddDistortions->isChecked(); + parameters.m_SignalGen.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()); ItkDoubleImgType::Pointer itkImg = ItkDoubleImgType::New(); CastToItkImage< ItkDoubleImgType >(img, itkImg); - if (parameters.m_ImageRegion.GetSize(0)==itkImg->GetLargestPossibleRegion().GetSize(0) && - parameters.m_ImageRegion.GetSize(1)==itkImg->GetLargestPossibleRegion().GetSize(1) && - parameters.m_ImageRegion.GetSize(2)==itkImg->GetLargestPossibleRegion().GetSize(2)) + if (parameters.m_SignalGen.m_ImageRegion.GetSize(0)==itkImg->GetLargestPossibleRegion().GetSize(0) && + parameters.m_SignalGen.m_ImageRegion.GetSize(1)==itkImg->GetLargestPossibleRegion().GetSize(1) && + parameters.m_SignalGen.m_ImageRegion.GetSize(2)==itkImg->GetLargestPossibleRegion().GetSize(2)) { - parameters.m_SimulateKspaceAcquisition = true; + parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; itk::ImageDuplicator::Pointer duplicator = itk::ImageDuplicator::New(); duplicator->SetInputImage(itkImg); duplicator->Update(); - parameters.m_FrequencyMap = duplicator->GetOutput(); - parameters.m_ArtifactModelString += "_DISTORTED"; - parameters.m_ResultNode->AddProperty("Fiberfox.Distortions", BoolProperty::New(true)); + parameters.m_SignalGen.m_FrequencyMap = duplicator->GetOutput(); + parameters.m_Misc.m_ArtifactModelString += "_DISTORTED"; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Distortions", BoolProperty::New(true)); } } - parameters.m_EddyStrength = 0; - parameters.m_DoAddEddyCurrents = m_Controls->m_AddEddy->isChecked(); + parameters.m_SignalGen.m_EddyStrength = 0; + parameters.m_SignalGen.m_DoAddEddyCurrents = m_Controls->m_AddEddy->isChecked(); if (m_Controls->m_AddEddy->isChecked()) { - parameters.m_EddyStrength = m_Controls->m_EddyGradientStrength->value(); - parameters.m_ArtifactModelString += "_EDDY"; - parameters.m_ResultNode->AddProperty("Fiberfox.Eddy-strength", DoubleProperty::New(parameters.m_EddyStrength)); + parameters.m_SignalGen.m_EddyStrength = m_Controls->m_EddyGradientStrength->value(); + parameters.m_Misc.m_ArtifactModelString += "_EDDY"; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Eddy-strength", DoubleProperty::New(parameters.m_SignalGen.m_EddyStrength)); } // Motion - parameters.m_DoAddMotion = m_Controls->m_AddMotion->isChecked(); - parameters.m_DoRandomizeMotion = m_Controls->m_RandomMotion->isChecked(); - parameters.m_Translation[0] = m_Controls->m_MaxTranslationBoxX->value(); - parameters.m_Translation[1] = m_Controls->m_MaxTranslationBoxY->value(); - parameters.m_Translation[2] = m_Controls->m_MaxTranslationBoxZ->value(); - parameters.m_Rotation[0] = m_Controls->m_MaxRotationBoxX->value(); - parameters.m_Rotation[1] = m_Controls->m_MaxRotationBoxY->value(); - parameters.m_Rotation[2] = m_Controls->m_MaxRotationBoxZ->value(); + parameters.m_SignalGen.m_DoAddMotion = m_Controls->m_AddMotion->isChecked(); + parameters.m_SignalGen.m_DoRandomizeMotion = m_Controls->m_RandomMotion->isChecked(); + parameters.m_SignalGen.m_Translation[0] = m_Controls->m_MaxTranslationBoxX->value(); + parameters.m_SignalGen.m_Translation[1] = m_Controls->m_MaxTranslationBoxY->value(); + parameters.m_SignalGen.m_Translation[2] = m_Controls->m_MaxTranslationBoxZ->value(); + parameters.m_SignalGen.m_Rotation[0] = m_Controls->m_MaxRotationBoxX->value(); + parameters.m_SignalGen.m_Rotation[1] = m_Controls->m_MaxRotationBoxY->value(); + parameters.m_SignalGen.m_Rotation[2] = m_Controls->m_MaxRotationBoxZ->value(); if ( m_Controls->m_AddMotion->isChecked() && m_SelectedBundles.size()>0 ) { - parameters.m_ArtifactModelString += "_MOTION"; - parameters.m_ResultNode->AddProperty("Fiberfox.Motion.Random", BoolProperty::New(parameters.m_DoRandomizeMotion)); - parameters.m_ResultNode->AddProperty("Fiberfox.Motion.Translation-x", DoubleProperty::New(parameters.m_Translation[0])); - parameters.m_ResultNode->AddProperty("Fiberfox.Motion.Translation-y", DoubleProperty::New(parameters.m_Translation[1])); - parameters.m_ResultNode->AddProperty("Fiberfox.Motion.Translation-z", DoubleProperty::New(parameters.m_Translation[2])); - parameters.m_ResultNode->AddProperty("Fiberfox.Motion.Rotation-x", DoubleProperty::New(parameters.m_Rotation[0])); - parameters.m_ResultNode->AddProperty("Fiberfox.Motion.Rotation-y", DoubleProperty::New(parameters.m_Rotation[1])); - parameters.m_ResultNode->AddProperty("Fiberfox.Motion.Rotation-z", DoubleProperty::New(parameters.m_Rotation[2])); + parameters.m_Misc.m_ArtifactModelString += "_MOTION"; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Random", BoolProperty::New(parameters.m_SignalGen.m_DoRandomizeMotion)); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Translation-x", DoubleProperty::New(parameters.m_SignalGen.m_Translation[0])); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Translation-y", DoubleProperty::New(parameters.m_SignalGen.m_Translation[1])); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Translation-z", DoubleProperty::New(parameters.m_SignalGen.m_Translation[2])); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Rotation-x", DoubleProperty::New(parameters.m_SignalGen.m_Rotation[0])); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Rotation-y", DoubleProperty::New(parameters.m_SignalGen.m_Rotation[1])); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Rotation-z", DoubleProperty::New(parameters.m_SignalGen.m_Rotation[2])); } // other imaging parameters - parameters.m_tLine = m_Controls->m_LineReadoutTimeBox->value(); - parameters.m_tInhom = m_Controls->m_T2starBox->value(); - parameters.m_tEcho = m_Controls->m_TEbox->value(); - parameters.m_DoDisablePartialVolume = m_Controls->m_EnforcePureFiberVoxelsBox->isChecked(); - parameters.m_AxonRadius = m_Controls->m_FiberRadius->value(); - parameters.m_SignalScale = m_Controls->m_SignalScaleBox->value(); + parameters.m_SignalGen.m_tLine = m_Controls->m_LineReadoutTimeBox->value(); + parameters.m_SignalGen.m_tInhom = m_Controls->m_T2starBox->value(); + parameters.m_SignalGen.m_tEcho = m_Controls->m_TEbox->value(); + parameters.m_SignalGen.m_DoDisablePartialVolume = m_Controls->m_EnforcePureFiberVoxelsBox->isChecked(); + parameters.m_SignalGen.m_AxonRadius = m_Controls->m_FiberRadius->value(); + parameters.m_SignalGen.m_SignalScale = m_Controls->m_SignalScaleBox->value(); // adjust echo time if needed - if ( parameters.m_tEcho < parameters.m_ImageRegion.GetSize(1)*parameters.m_tLine ) + if ( parameters.m_SignalGen.m_tEcho < parameters.m_SignalGen.m_ImageRegion.GetSize(1)*parameters.m_SignalGen.m_tLine ) { - this->m_Controls->m_TEbox->setValue( parameters.m_ImageRegion.GetSize(1)*parameters.m_tLine ); - parameters.m_tEcho = m_Controls->m_TEbox->value(); - QMessageBox::information( NULL, "Warning", "Echo time is too short! Time not sufficient to read slice. Automaticall adjusted to "+QString::number(parameters.m_tEcho)+" ms"); + this->m_Controls->m_TEbox->setValue( parameters.m_SignalGen.m_ImageRegion.GetSize(1)*parameters.m_SignalGen.m_tLine ); + parameters.m_SignalGen.m_tEcho = m_Controls->m_TEbox->value(); + QMessageBox::information( NULL, "Warning", "Echo time is too short! Time not sufficient to read slice. Automaticall adjusted to "+QString::number(parameters.m_SignalGen.m_tEcho)+" ms"); } // Noise - parameters.m_DoAddNoise = m_Controls->m_AddNoise->isChecked(); + parameters.m_SignalGen.m_DoAddNoise = m_Controls->m_AddNoise->isChecked(); if (m_Controls->m_AddNoise->isChecked()) { double noiseVariance = m_Controls->m_NoiseLevel->value(); { switch (m_Controls->m_NoiseDistributionBox->currentIndex()) { case 0: { parameters.m_NoiseModel = new mitk::RicianNoiseModel(); - parameters.m_ArtifactModelString += "_RICIAN-"; - parameters.m_ResultNode->AddProperty("Fiberfox.Noise-Distribution", StringProperty::New("Rician")); + parameters.m_Misc.m_ArtifactModelString += "_RICIAN-"; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Noise-Distribution", StringProperty::New("Rician")); break; } case 1: { parameters.m_NoiseModel = new mitk::ChiSquareNoiseModel(); - parameters.m_ArtifactModelString += "_CHISQUARED-"; - parameters.m_ResultNode->AddProperty("Fiberfox.Noise-Distribution", StringProperty::New("Chi-squared")); + parameters.m_Misc.m_ArtifactModelString += "_CHISQUARED-"; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Noise-Distribution", StringProperty::New("Chi-squared")); break; } default: { parameters.m_NoiseModel = new mitk::RicianNoiseModel(); - parameters.m_ArtifactModelString += "_RICIAN-"; - parameters.m_ResultNode->AddProperty("Fiberfox.Noise-Distribution", StringProperty::New("Rician")); + parameters.m_Misc.m_ArtifactModelString += "_RICIAN-"; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Noise-Distribution", StringProperty::New("Rician")); } } } parameters.m_NoiseModel->SetNoiseVariance(noiseVariance); - parameters.m_ArtifactModelString += QString::number(noiseVariance).toStdString(); - parameters.m_ResultNode->AddProperty("Fiberfox.Noise-Variance", DoubleProperty::New(noiseVariance)); + parameters.m_Misc.m_ArtifactModelString += QString::number(noiseVariance).toStdString(); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Noise-Variance", DoubleProperty::New(noiseVariance)); } // adjusting line readout time to the adapted image size needed for the DFT - unsigned int y = parameters.m_ImageRegion.GetSize(1); + unsigned int y = parameters.m_SignalGen.m_ImageRegion.GetSize(1); y += y%2; - if ( y>parameters.m_ImageRegion.GetSize(1) ) - parameters.m_tLine *= (double)parameters.m_ImageRegion.GetSize(1)/y; + if ( y>parameters.m_SignalGen.m_ImageRegion.GetSize(1) ) + parameters.m_SignalGen.m_tLine *= (double)parameters.m_SignalGen.m_ImageRegion.GetSize(1)/y; // signal models m_PrototypeModel1.Clear(); m_PrototypeModel3.Clear(); m_PrototypeModel4.Clear(); // compartment 1 switch (m_Controls->m_Compartment1Box->currentIndex()) { case 0: - m_StickModel1.SetGradientList(parameters.GetGradientDirections()); - m_StickModel1.SetBvalue(parameters.m_Bvalue); - m_StickModel1.SetDiffusivity(m_Controls->m_StickWidget1->GetD()); - m_StickModel1.SetT2(m_Controls->m_StickWidget1->GetT2()); - parameters.m_FiberModelList.push_back(&m_StickModel1); - parameters.m_SignalModelString += "Stick"; - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Stick") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment1.D", DoubleProperty::New(m_Controls->m_StickWidget1->GetD()) ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment1.T2", DoubleProperty::New(m_StickModel1.GetT2()) ); + mitk::StickModel* model = new mitk::StickModel(); + model->SetGradientList(parameters.GetGradientDirections()); + model->SetBvalue(parameters.m_SignalGen.m_Bvalue); + model->SetDiffusivity(m_Controls->m_StickWidget1->GetD()); + model->SetT2(m_Controls->m_StickWidget1->GetT2()); + parameters.m_FiberModelList.push_back(model); + parameters.m_Misc.m_SignalModelString += "Stick"; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Stick") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D", DoubleProperty::New(m_Controls->m_StickWidget1->GetD()) ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.T2", DoubleProperty::New(m_StickModel1.GetT2()) ); break; case 1: m_ZeppelinModel1.SetGradientList(parameters.GetGradientDirections()); - m_ZeppelinModel1.SetBvalue(parameters.m_Bvalue); + m_ZeppelinModel1.SetBvalue(parameters.m_SignalGen.m_Bvalue); m_ZeppelinModel1.SetDiffusivity1(m_Controls->m_ZeppelinWidget1->GetD1()); m_ZeppelinModel1.SetDiffusivity2(m_Controls->m_ZeppelinWidget1->GetD2()); m_ZeppelinModel1.SetDiffusivity3(m_Controls->m_ZeppelinWidget1->GetD2()); m_ZeppelinModel1.SetT2(m_Controls->m_ZeppelinWidget1->GetT2()); parameters.m_FiberModelList.push_back(&m_ZeppelinModel1); - parameters.m_SignalModelString += "Zeppelin"; - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Zeppelin") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment1.D1", DoubleProperty::New(m_Controls->m_ZeppelinWidget1->GetD1()) ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment1.D2", DoubleProperty::New(m_Controls->m_ZeppelinWidget1->GetD2()) ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment1.T2", DoubleProperty::New(m_ZeppelinModel1.GetT2()) ); + parameters.m_Misc.m_SignalModelString += "Zeppelin"; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Zeppelin") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D1", DoubleProperty::New(m_Controls->m_ZeppelinWidget1->GetD1()) ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D2", DoubleProperty::New(m_Controls->m_ZeppelinWidget1->GetD2()) ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.T2", DoubleProperty::New(m_ZeppelinModel1.GetT2()) ); break; case 2: m_TensorModel1.SetGradientList(parameters.GetGradientDirections()); - m_TensorModel1.SetBvalue(parameters.m_Bvalue); + m_TensorModel1.SetBvalue(parameters.m_SignalGen.m_Bvalue); m_TensorModel1.SetDiffusivity1(m_Controls->m_TensorWidget1->GetD1()); m_TensorModel1.SetDiffusivity2(m_Controls->m_TensorWidget1->GetD2()); m_TensorModel1.SetDiffusivity3(m_Controls->m_TensorWidget1->GetD3()); m_TensorModel1.SetT2(m_Controls->m_TensorWidget1->GetT2()); parameters.m_FiberModelList.push_back(&m_TensorModel1); - parameters.m_SignalModelString += "Tensor"; - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Tensor") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment1.D1", DoubleProperty::New(m_Controls->m_TensorWidget1->GetD1()) ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment1.D2", DoubleProperty::New(m_Controls->m_TensorWidget1->GetD2()) ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment1.D3", DoubleProperty::New(m_Controls->m_TensorWidget1->GetD3()) ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment1.T2", DoubleProperty::New(m_ZeppelinModel1.GetT2()) ); + parameters.m_Misc.m_SignalModelString += "Tensor"; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Tensor") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D1", DoubleProperty::New(m_Controls->m_TensorWidget1->GetD1()) ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D2", DoubleProperty::New(m_Controls->m_TensorWidget1->GetD2()) ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D3", DoubleProperty::New(m_Controls->m_TensorWidget1->GetD3()) ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.T2", DoubleProperty::New(m_ZeppelinModel1.GetT2()) ); break; case 3: - parameters.m_SimulateKspaceAcquisition = false; + parameters.m_SignalGen.m_SimulateKspaceAcquisition = false; m_PrototypeModel1.SetGradientList(parameters.GetGradientDirections()); m_PrototypeModel1.SetMaxNumKernels(m_Controls->m_PrototypeWidget1->GetNumberOfSamples()); m_PrototypeModel1.SetFaRange(m_Controls->m_PrototypeWidget1->GetMinFa(), m_Controls->m_PrototypeWidget1->GetMaxFa()); m_PrototypeModel1.SetAdcRange(m_Controls->m_PrototypeWidget1->GetMinAdc(), m_Controls->m_PrototypeWidget1->GetMaxAdc()); parameters.m_FiberModelList.push_back(&m_PrototypeModel1); - parameters.m_SignalModelString += "Prototype"; - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Prototype") ); + parameters.m_Misc.m_SignalModelString += "Prototype"; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Prototype") ); break; } // compartment 2 switch (m_Controls->m_Compartment2Box->currentIndex()) { case 0: break; case 1: m_StickModel2.SetGradientList(parameters.GetGradientDirections()); - m_StickModel2.SetBvalue(parameters.m_Bvalue); + m_StickModel2.SetBvalue(parameters.m_SignalGen.m_Bvalue); m_StickModel2.SetDiffusivity(m_Controls->m_StickWidget2->GetD()); m_StickModel2.SetT2(m_Controls->m_StickWidget2->GetT2()); parameters.m_FiberModelList.push_back(&m_StickModel2); - parameters.m_SignalModelString += "Stick"; - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment2.Description", StringProperty::New("Inter-axonal compartment") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment2.Model", StringProperty::New("Stick") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment2.D", DoubleProperty::New(m_Controls->m_StickWidget2->GetD()) ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment2.T2", DoubleProperty::New(m_StickModel2.GetT2()) ); + parameters.m_Misc.m_SignalModelString += "Stick"; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Description", StringProperty::New("Inter-axonal compartment") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Model", StringProperty::New("Stick") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D", DoubleProperty::New(m_Controls->m_StickWidget2->GetD()) ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.T2", DoubleProperty::New(m_StickModel2.GetT2()) ); break; case 2: m_ZeppelinModel2.SetGradientList(parameters.GetGradientDirections()); - m_ZeppelinModel2.SetBvalue(parameters.m_Bvalue); + m_ZeppelinModel2.SetBvalue(parameters.m_SignalGen.m_Bvalue); m_ZeppelinModel2.SetDiffusivity1(m_Controls->m_ZeppelinWidget2->GetD1()); m_ZeppelinModel2.SetDiffusivity2(m_Controls->m_ZeppelinWidget2->GetD2()); m_ZeppelinModel2.SetDiffusivity3(m_Controls->m_ZeppelinWidget2->GetD2()); m_ZeppelinModel2.SetT2(m_Controls->m_ZeppelinWidget2->GetT2()); parameters.m_FiberModelList.push_back(&m_ZeppelinModel2); - parameters.m_SignalModelString += "Zeppelin"; - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment2.Description", StringProperty::New("Inter-axonal compartment") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment2.Model", StringProperty::New("Zeppelin") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment2.D1", DoubleProperty::New(m_Controls->m_ZeppelinWidget2->GetD1()) ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment2.D2", DoubleProperty::New(m_Controls->m_ZeppelinWidget2->GetD2()) ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment2.T2", DoubleProperty::New(m_ZeppelinModel2.GetT2()) ); + parameters.m_Misc.m_SignalModelString += "Zeppelin"; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Description", StringProperty::New("Inter-axonal compartment") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Model", StringProperty::New("Zeppelin") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D1", DoubleProperty::New(m_Controls->m_ZeppelinWidget2->GetD1()) ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D2", DoubleProperty::New(m_Controls->m_ZeppelinWidget2->GetD2()) ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.T2", DoubleProperty::New(m_ZeppelinModel2.GetT2()) ); break; case 3: m_TensorModel2.SetGradientList(parameters.GetGradientDirections()); - m_TensorModel2.SetBvalue(parameters.m_Bvalue); + m_TensorModel2.SetBvalue(parameters.m_SignalGen.m_Bvalue); m_TensorModel2.SetDiffusivity1(m_Controls->m_TensorWidget2->GetD1()); m_TensorModel2.SetDiffusivity2(m_Controls->m_TensorWidget2->GetD2()); m_TensorModel2.SetDiffusivity3(m_Controls->m_TensorWidget2->GetD3()); m_TensorModel2.SetT2(m_Controls->m_TensorWidget2->GetT2()); parameters.m_FiberModelList.push_back(&m_TensorModel2); - parameters.m_SignalModelString += "Tensor"; - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment2.Description", StringProperty::New("Inter-axonal compartment") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment2.Model", StringProperty::New("Tensor") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment2.D1", DoubleProperty::New(m_Controls->m_TensorWidget2->GetD1()) ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment2.D2", DoubleProperty::New(m_Controls->m_TensorWidget2->GetD2()) ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment2.D3", DoubleProperty::New(m_Controls->m_TensorWidget2->GetD3()) ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment2.T2", DoubleProperty::New(m_ZeppelinModel2.GetT2()) ); + parameters.m_Misc.m_SignalModelString += "Tensor"; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Description", StringProperty::New("Inter-axonal compartment") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Model", StringProperty::New("Tensor") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D1", DoubleProperty::New(m_Controls->m_TensorWidget2->GetD1()) ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D2", DoubleProperty::New(m_Controls->m_TensorWidget2->GetD2()) ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D3", DoubleProperty::New(m_Controls->m_TensorWidget2->GetD3()) ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.T2", DoubleProperty::New(m_ZeppelinModel2.GetT2()) ); break; } // compartment 3 switch (m_Controls->m_Compartment3Box->currentIndex()) { case 0: m_BallModel1.SetGradientList(parameters.GetGradientDirections()); - m_BallModel1.SetBvalue(parameters.m_Bvalue); + m_BallModel1.SetBvalue(parameters.m_SignalGen.m_Bvalue); m_BallModel1.SetDiffusivity(m_Controls->m_BallWidget1->GetD()); m_BallModel1.SetT2(m_Controls->m_BallWidget1->GetT2()); parameters.m_NonFiberModelList.push_back(&m_BallModel1); - parameters.m_SignalModelString += "Ball"; - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Ball") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment3.D", DoubleProperty::New(m_Controls->m_BallWidget1->GetD()) ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment3.T2", DoubleProperty::New(m_BallModel1.GetT2()) ); + parameters.m_Misc.m_SignalModelString += "Ball"; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Ball") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.D", DoubleProperty::New(m_Controls->m_BallWidget1->GetD()) ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.T2", DoubleProperty::New(m_BallModel1.GetT2()) ); break; case 1: m_AstrosticksModel1.SetGradientList(parameters.GetGradientDirections()); - m_AstrosticksModel1.SetBvalue(parameters.m_Bvalue); + m_AstrosticksModel1.SetBvalue(parameters.m_SignalGen.m_Bvalue); m_AstrosticksModel1.SetDiffusivity(m_Controls->m_AstrosticksWidget1->GetD()); m_AstrosticksModel1.SetT2(m_Controls->m_AstrosticksWidget1->GetT2()); m_AstrosticksModel1.SetRandomizeSticks(m_Controls->m_AstrosticksWidget1->GetRandomizeSticks()); parameters.m_NonFiberModelList.push_back(&m_AstrosticksModel1); - parameters.m_SignalModelString += "Astrosticks"; - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Astrosticks") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment3.D", DoubleProperty::New(m_Controls->m_AstrosticksWidget1->GetD()) ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment3.T2", DoubleProperty::New(m_AstrosticksModel1.GetT2()) ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment3.RandomSticks", BoolProperty::New(m_Controls->m_AstrosticksWidget1->GetRandomizeSticks()) ); + parameters.m_Misc.m_SignalModelString += "Astrosticks"; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Astrosticks") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.D", DoubleProperty::New(m_Controls->m_AstrosticksWidget1->GetD()) ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.T2", DoubleProperty::New(m_AstrosticksModel1.GetT2()) ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.RandomSticks", BoolProperty::New(m_Controls->m_AstrosticksWidget1->GetRandomizeSticks()) ); break; case 2: m_DotModel1.SetGradientList(parameters.GetGradientDirections()); m_DotModel1.SetT2(m_Controls->m_DotWidget1->GetT2()); parameters.m_NonFiberModelList.push_back(&m_DotModel1); - parameters.m_SignalModelString += "Dot"; - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Dot") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment3.T2", DoubleProperty::New(m_DotModel1.GetT2()) ); + parameters.m_Misc.m_SignalModelString += "Dot"; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Dot") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.T2", DoubleProperty::New(m_DotModel1.GetT2()) ); break; case 3: - parameters.m_SimulateKspaceAcquisition = false; + parameters.m_SignalGen.m_SimulateKspaceAcquisition = false; m_PrototypeModel3.SetGradientList(parameters.GetGradientDirections()); m_PrototypeModel3.SetMaxNumKernels(m_Controls->m_PrototypeWidget3->GetNumberOfSamples()); m_PrototypeModel3.SetFaRange(m_Controls->m_PrototypeWidget3->GetMinFa(), m_Controls->m_PrototypeWidget3->GetMaxFa()); m_PrototypeModel3.SetAdcRange(m_Controls->m_PrototypeWidget3->GetMinAdc(), m_Controls->m_PrototypeWidget3->GetMaxAdc()); parameters.m_NonFiberModelList.push_back(&m_PrototypeModel3); - parameters.m_SignalModelString += "Prototype"; - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Prototype") ); + parameters.m_Misc.m_SignalModelString += "Prototype"; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Prototype") ); break; } // compartment 4 ItkDoubleImgType::Pointer comp4VolumeImage = NULL; ItkDoubleImgType::Pointer comp3VolumeImage = NULL; if (m_Controls->m_Compartment4Box->currentIndex()>0) { mitk::DataNode::Pointer volumeNode = m_Controls->m_Comp4VolumeFraction->GetSelectedNode(); if (volumeNode.IsNull()) { QMessageBox::information( NULL, "Information", "No volume fraction image selected! Second extra-axonal compartment has been disabled for this simultation."); MITK_WARN << "No volume fraction image selected! Second extra-axonal compartment has been disabled."; } else { MITK_INFO << "Rescaling volume fraction image..."; comp4VolumeImage = ItkDoubleImgType::New(); mitk::Image* img = dynamic_cast(volumeNode->GetData()); CastToItkImage< ItkDoubleImgType >(img, comp4VolumeImage); double max = itk::NumericTraits::min(); double min = itk::NumericTraits::max(); itk::ImageRegionIterator< ItkDoubleImgType > it(comp4VolumeImage, comp4VolumeImage->GetLargestPossibleRegion()); while(!it.IsAtEnd()) { - if (parameters.m_MaskImage.IsNotNull() && parameters.m_MaskImage->GetPixel(it.GetIndex())<=0) + if (parameters.m_SignalGen.m_MaskImage.IsNotNull() && 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()::Pointer scaler = itk::ShiftScaleImageFilter< ItkDoubleImgType, ItkDoubleImgType >::New(); scaler->SetInput(comp4VolumeImage); scaler->SetShift(-min); scaler->SetScale(1.0/(max-min)); scaler->Update(); comp4VolumeImage = scaler->GetOutput(); // itk::ImageFileWriter< ItkDoubleImgType >::Pointer wr = itk::ImageFileWriter< ItkDoubleImgType >::New(); // wr->SetInput(comp4VolumeImage); // wr->SetFileName("/local/comp4.nrrd"); // wr->Update(); // if (max>1 || min<0) // are volume fractions between 0 and 1? // { // itk::RescaleIntensityImageFilter::Pointer rescaler = itk::RescaleIntensityImageFilter::New(); // rescaler->SetInput(0, comp4VolumeImage); // rescaler->SetOutputMaximum(1); // rescaler->SetOutputMinimum(0); // rescaler->Update(); // comp4VolumeImage = rescaler->GetOutput(); // } itk::InvertIntensityImageFilter< ItkDoubleImgType, ItkDoubleImgType >::Pointer inverter = itk::InvertIntensityImageFilter< ItkDoubleImgType, ItkDoubleImgType >::New(); inverter->SetMaximum(1.0); inverter->SetInput(comp4VolumeImage); inverter->Update(); comp3VolumeImage = inverter->GetOutput(); } } if (comp4VolumeImage.IsNotNull()) switch (m_Controls->m_Compartment4Box->currentIndex()) { case 0: break; case 1: { m_BallModel2.SetGradientList(parameters.GetGradientDirections()); - m_BallModel2.SetBvalue(parameters.m_Bvalue); + m_BallModel2.SetBvalue(parameters.m_SignalGen.m_Bvalue); m_BallModel2.SetDiffusivity(m_Controls->m_BallWidget2->GetD()); m_BallModel2.SetT2(m_Controls->m_BallWidget2->GetT2()); m_BallModel2.SetVolumeFractionImage(comp4VolumeImage); parameters.m_NonFiberModelList.back()->SetVolumeFractionImage(comp3VolumeImage); parameters.m_NonFiberModelList.push_back(&m_BallModel2); - parameters.m_SignalModelString += "Ball"; - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Ball") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment4.D", DoubleProperty::New(m_Controls->m_BallWidget2->GetD()) ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment4.T2", DoubleProperty::New(m_BallModel2.GetT2()) ); + parameters.m_Misc.m_SignalModelString += "Ball"; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Ball") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.D", DoubleProperty::New(m_Controls->m_BallWidget2->GetD()) ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.T2", DoubleProperty::New(m_BallModel2.GetT2()) ); break; } case 2: { m_AstrosticksModel2.SetGradientList(parameters.GetGradientDirections()); - m_AstrosticksModel2.SetBvalue(parameters.m_Bvalue); + m_AstrosticksModel2.SetBvalue(parameters.m_SignalGen.m_Bvalue); m_AstrosticksModel2.SetDiffusivity(m_Controls->m_AstrosticksWidget2->GetD()); m_AstrosticksModel2.SetT2(m_Controls->m_AstrosticksWidget2->GetT2()); m_AstrosticksModel2.SetRandomizeSticks(m_Controls->m_AstrosticksWidget2->GetRandomizeSticks()); parameters.m_NonFiberModelList.back()->SetVolumeFractionImage(comp3VolumeImage); m_AstrosticksModel2.SetVolumeFractionImage(comp4VolumeImage); parameters.m_NonFiberModelList.push_back(&m_AstrosticksModel2); - parameters.m_SignalModelString += "Astrosticks"; - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Astrosticks") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment4.D", DoubleProperty::New(m_Controls->m_AstrosticksWidget2->GetD()) ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment4.T2", DoubleProperty::New(m_AstrosticksModel2.GetT2()) ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment4.RandomSticks", BoolProperty::New(m_Controls->m_AstrosticksWidget2->GetRandomizeSticks()) ); + parameters.m_Misc.m_SignalModelString += "Astrosticks"; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Astrosticks") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.D", DoubleProperty::New(m_Controls->m_AstrosticksWidget2->GetD()) ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.T2", DoubleProperty::New(m_AstrosticksModel2.GetT2()) ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.RandomSticks", BoolProperty::New(m_Controls->m_AstrosticksWidget2->GetRandomizeSticks()) ); break; } case 3: { m_DotModel2.SetGradientList(parameters.GetGradientDirections()); m_DotModel2.SetT2(m_Controls->m_DotWidget2->GetT2()); m_DotModel2.SetVolumeFractionImage(comp4VolumeImage); parameters.m_NonFiberModelList.back()->SetVolumeFractionImage(comp3VolumeImage); parameters.m_NonFiberModelList.push_back(&m_DotModel2); - parameters.m_SignalModelString += "Dot"; - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Dot") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment4.T2", DoubleProperty::New(m_DotModel2.GetT2()) ); + parameters.m_Misc.m_SignalModelString += "Dot"; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Dot") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.T2", DoubleProperty::New(m_DotModel2.GetT2()) ); break; } case 4: { - parameters.m_SimulateKspaceAcquisition = false; + parameters.m_SignalGen.m_SimulateKspaceAcquisition = false; m_PrototypeModel4.SetGradientList(parameters.GetGradientDirections()); m_PrototypeModel4.SetMaxNumKernels(m_Controls->m_PrototypeWidget4->GetNumberOfSamples()); m_PrototypeModel4.SetFaRange(m_Controls->m_PrototypeWidget4->GetMinFa(), m_Controls->m_PrototypeWidget4->GetMaxFa()); m_PrototypeModel4.SetAdcRange(m_Controls->m_PrototypeWidget4->GetMinAdc(), m_Controls->m_PrototypeWidget4->GetMaxAdc()); m_PrototypeModel4.SetVolumeFractionImage(comp4VolumeImage); parameters.m_NonFiberModelList.back()->SetVolumeFractionImage(comp3VolumeImage); parameters.m_NonFiberModelList.push_back(&m_PrototypeModel4); - parameters.m_SignalModelString += "Prototype"; - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); - parameters.m_ResultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Prototype") ); + parameters.m_Misc.m_SignalModelString += "Prototype"; + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Prototype") ); break; } } - parameters.m_FiberSeparationThreshold = m_Controls->m_SeparationAngleBox->value(); + parameters.m_SignalGen.m_FiberSeparationThreshold = m_Controls->m_SeparationAngleBox->value(); switch (m_Controls->m_DiffusionDirectionBox->currentIndex()) { case 0: - parameters.m_DiffusionDirectionMode = FiberfoxParameters::FIBER_TANGENT_DIRECTIONS; + parameters.m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::FIBER_TANGENT_DIRECTIONS; break; case 1: - parameters.m_DiffusionDirectionMode = FiberfoxParameters::MAIN_FIBER_DIRECTIONS; + parameters.m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::MAIN_FIBER_DIRECTIONS; break; case 2: - parameters.m_DiffusionDirectionMode = FiberfoxParameters::RANDOM_DIRECTIONS; - parameters.m_DoAddMotion = false; - parameters.m_DoAddGibbsRinging = false; - parameters.m_KspaceLineOffset = 0.0; - parameters.m_FrequencyMap = NULL; - parameters.m_CroppingFactor = 1.0; - parameters.m_EddyStrength = 0; + parameters.m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::RANDOM_DIRECTIONS; + parameters.m_SignalGen.m_DoAddMotion = false; + parameters.m_SignalGen.m_DoAddGibbsRinging = false; + parameters.m_SignalGen.m_KspaceLineOffset = 0.0; + parameters.m_SignalGen.m_FrequencyMap = NULL; + parameters.m_SignalGen.m_CroppingFactor = 1.0; + parameters.m_SignalGen.m_EddyStrength = 0; break; default: - parameters.m_DiffusionDirectionMode = FiberfoxParameters::FIBER_TANGENT_DIRECTIONS; - } - - parameters.m_ResultNode->AddProperty("Fiberfox.SignalScale", IntProperty::New(parameters.m_SignalScale)); - parameters.m_ResultNode->AddProperty("Fiberfox.FiberRadius", IntProperty::New(parameters.m_AxonRadius)); - parameters.m_ResultNode->AddProperty("Fiberfox.Tinhom", DoubleProperty::New(parameters.m_tInhom)); - parameters.m_ResultNode->AddProperty("Fiberfox.Tline", DoubleProperty::New(parameters.m_tLine)); - parameters.m_ResultNode->AddProperty("Fiberfox.TE", DoubleProperty::New(parameters.m_tEcho)); - parameters.m_ResultNode->AddProperty("Fiberfox.b-value", DoubleProperty::New(parameters.m_Bvalue)); - parameters.m_ResultNode->AddProperty("Fiberfox.NoPartialVolume", BoolProperty::New(parameters.m_DoDisablePartialVolume)); - parameters.m_ResultNode->AddProperty("Fiberfox.Relaxation", BoolProperty::New(parameters.m_DoSimulateRelaxation)); - parameters.m_ResultNode->AddProperty("binary", BoolProperty::New(false)); - - parameters.m_FiberGenerationParameters.m_RealTimeFibers = m_Controls->m_RealTimeFibers->isChecked(); - parameters.m_FiberGenerationParameters.m_AdvancedOptions = m_Controls->m_AdvancedOptionsBox->isChecked(); - parameters.m_FiberGenerationParameters.m_Distribution = m_Controls->m_DistributionBox->currentIndex(); - parameters.m_FiberGenerationParameters.m_Variance = m_Controls->m_VarianceBox->value(); - parameters.m_FiberGenerationParameters.m_FiberDensity = m_Controls->m_FiberDensityBox->value(); - parameters.m_FiberGenerationParameters.m_IncludeFiducials = m_Controls->m_IncludeFiducials->isChecked(); - parameters.m_FiberGenerationParameters.m_ConstantRadius = m_Controls->m_ConstantRadiusBox->isChecked(); - parameters.m_FiberGenerationParameters.m_Sampling = m_Controls->m_FiberSamplingBox->value(); - parameters.m_FiberGenerationParameters.m_Tension = m_Controls->m_TensionBox->value(); - parameters.m_FiberGenerationParameters.m_Continuity = m_Controls->m_ContinuityBox->value(); - parameters.m_FiberGenerationParameters.m_Bias = m_Controls->m_BiasBox->value(); - parameters.m_FiberGenerationParameters.m_Rotation[0] = m_Controls->m_XrotBox->value(); - parameters.m_FiberGenerationParameters.m_Rotation[1] = m_Controls->m_YrotBox->value(); - parameters.m_FiberGenerationParameters.m_Rotation[2] = m_Controls->m_ZrotBox->value(); - parameters.m_FiberGenerationParameters.m_Translation[0] = m_Controls->m_XtransBox->value(); - parameters.m_FiberGenerationParameters.m_Translation[1] = m_Controls->m_YtransBox->value(); - parameters.m_FiberGenerationParameters.m_Translation[2] = m_Controls->m_ZtransBox->value(); - parameters.m_FiberGenerationParameters.m_Scale[0] = m_Controls->m_XscaleBox->value(); - parameters.m_FiberGenerationParameters.m_Scale[1] = m_Controls->m_YscaleBox->value(); - parameters.m_FiberGenerationParameters.m_Scale[2] = m_Controls->m_ZscaleBox->value(); + parameters.m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::FIBER_TANGENT_DIRECTIONS; + } + + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.SignalScale", IntProperty::New(parameters.m_SignalGen.m_SignalScale)); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.FiberRadius", IntProperty::New(parameters.m_SignalGen.m_AxonRadius)); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Tinhom", DoubleProperty::New(parameters.m_SignalGen.m_tInhom)); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Tline", DoubleProperty::New(parameters.m_SignalGen.m_tLine)); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.TE", DoubleProperty::New(parameters.m_SignalGen.m_tEcho)); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.b-value", DoubleProperty::New(parameters.m_SignalGen.m_Bvalue)); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.NoPartialVolume", BoolProperty::New(parameters.m_SignalGen.m_DoDisablePartialVolume)); + parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Relaxation", BoolProperty::New(parameters.m_SignalGen.m_DoSimulateRelaxation)); + parameters.m_Misc.m_ResultNode->AddProperty("binary", BoolProperty::New(false)); + + parameters.m_FiberGen.m_RealTimeFibers = m_Controls->m_RealTimeFibers->isChecked(); + parameters.m_FiberGen.m_AdvancedOptions = m_Controls->m_AdvancedOptionsBox->isChecked(); + parameters.m_FiberGen.m_Distribution = m_Controls->m_DistributionBox->currentIndex(); + parameters.m_FiberGen.m_Variance = m_Controls->m_VarianceBox->value(); + parameters.m_FiberGen.m_FiberDensity = m_Controls->m_FiberDensityBox->value(); + parameters.m_FiberGen.m_IncludeFiducials = m_Controls->m_IncludeFiducials->isChecked(); + parameters.m_FiberGen.m_ConstantRadius = m_Controls->m_ConstantRadiusBox->isChecked(); + parameters.m_FiberGen.m_Sampling = m_Controls->m_FiberSamplingBox->value(); + parameters.m_FiberGen.m_Tension = m_Controls->m_TensionBox->value(); + parameters.m_FiberGen.m_Continuity = m_Controls->m_ContinuityBox->value(); + parameters.m_FiberGen.m_Bias = m_Controls->m_BiasBox->value(); + parameters.m_FiberGen.m_Rotation[0] = m_Controls->m_XrotBox->value(); + parameters.m_FiberGen.m_Rotation[1] = m_Controls->m_YrotBox->value(); + parameters.m_FiberGen.m_Rotation[2] = m_Controls->m_ZrotBox->value(); + parameters.m_FiberGen.m_Translation[0] = m_Controls->m_XtransBox->value(); + parameters.m_FiberGen.m_Translation[1] = m_Controls->m_YtransBox->value(); + parameters.m_FiberGen.m_Translation[2] = m_Controls->m_ZtransBox->value(); + parameters.m_FiberGen.m_Scale[0] = m_Controls->m_XscaleBox->value(); + parameters.m_FiberGen.m_Scale[1] = m_Controls->m_YscaleBox->value(); + parameters.m_FiberGen.m_Scale[2] = m_Controls->m_ZscaleBox->value(); return parameters; } void QmitkFiberfoxView::SaveParameters() { FiberfoxParameters<> ffParamaters = UpdateImageParameters(); QString filename = QFileDialog::getSaveFileName( 0, tr("Save Parameters"), m_ParameterFile, tr("Fiberfox Parameters (*.ffp)") ); ffParamaters.SaveParameters(filename.toStdString()); m_ParameterFile = filename; } void QmitkFiberfoxView::LoadParameters() { QString filename = QFileDialog::getOpenFileName(0, tr("Load Parameters"), QString(itksys::SystemTools::GetFilenamePath(m_ParameterFile.toStdString()).c_str()), tr("Fiberfox Parameters (*.ffp)") ); if(filename.isEmpty() || filename.isNull()) return; m_ParameterFile = filename; - FiberfoxParameters<> ffParamaters; - ffParamaters.LoadParameters(filename.toStdString()); - - m_Controls->m_RealTimeFibers->setChecked(ffParamaters.m_FiberGenerationParameters.m_RealTimeFibers); - m_Controls->m_AdvancedOptionsBox->setChecked(ffParamaters.m_FiberGenerationParameters.m_AdvancedOptions); - m_Controls->m_DistributionBox->setCurrentIndex(ffParamaters.m_FiberGenerationParameters.m_Distribution); - m_Controls->m_VarianceBox->setValue(ffParamaters.m_FiberGenerationParameters.m_Variance); - m_Controls->m_FiberDensityBox->setValue(ffParamaters.m_FiberGenerationParameters.m_FiberDensity); - m_Controls->m_IncludeFiducials->setChecked(ffParamaters.m_FiberGenerationParameters.m_IncludeFiducials); - m_Controls->m_ConstantRadiusBox->setChecked(ffParamaters.m_FiberGenerationParameters.m_ConstantRadius); - m_Controls->m_FiberSamplingBox->setValue(ffParamaters.m_FiberGenerationParameters.m_Sampling); - m_Controls->m_TensionBox->setValue(ffParamaters.m_FiberGenerationParameters.m_Tension); - m_Controls->m_ContinuityBox->setValue(ffParamaters.m_FiberGenerationParameters.m_Continuity); - m_Controls->m_BiasBox->setValue(ffParamaters.m_FiberGenerationParameters.m_Bias); - m_Controls->m_XrotBox->setValue(ffParamaters.m_FiberGenerationParameters.m_Rotation[0]); - m_Controls->m_YrotBox->setValue(ffParamaters.m_FiberGenerationParameters.m_Rotation[1]); - m_Controls->m_ZrotBox->setValue(ffParamaters.m_FiberGenerationParameters.m_Rotation[2]); - m_Controls->m_XtransBox->setValue(ffParamaters.m_FiberGenerationParameters.m_Translation[0]); - m_Controls->m_YtransBox->setValue(ffParamaters.m_FiberGenerationParameters.m_Translation[1]); - m_Controls->m_ZtransBox->setValue(ffParamaters.m_FiberGenerationParameters.m_Translation[2]); - m_Controls->m_XscaleBox->setValue(ffParamaters.m_FiberGenerationParameters.m_Scale[0]); - m_Controls->m_YscaleBox->setValue(ffParamaters.m_FiberGenerationParameters.m_Scale[1]); - m_Controls->m_ZscaleBox->setValue(ffParamaters.m_FiberGenerationParameters.m_Scale[2]); + FiberfoxParameters<> parameters; + parameters.LoadParameters(filename.toStdString()); + + m_Controls->m_RealTimeFibers->setChecked(parameters.m_FiberGen.m_RealTimeFibers); + m_Controls->m_AdvancedOptionsBox->setChecked(parameters.m_FiberGen.m_AdvancedOptions); + m_Controls->m_DistributionBox->setCurrentIndex(parameters.m_FiberGen.m_Distribution); + m_Controls->m_VarianceBox->setValue(parameters.m_FiberGen.m_Variance); + m_Controls->m_FiberDensityBox->setValue(parameters.m_FiberGen.m_FiberDensity); + m_Controls->m_IncludeFiducials->setChecked(parameters.m_FiberGen.m_IncludeFiducials); + m_Controls->m_ConstantRadiusBox->setChecked(parameters.m_FiberGen.m_ConstantRadius); + m_Controls->m_FiberSamplingBox->setValue(parameters.m_FiberGen.m_Sampling); + m_Controls->m_TensionBox->setValue(parameters.m_FiberGen.m_Tension); + m_Controls->m_ContinuityBox->setValue(parameters.m_FiberGen.m_Continuity); + m_Controls->m_BiasBox->setValue(parameters.m_FiberGen.m_Bias); + m_Controls->m_XrotBox->setValue(parameters.m_FiberGen.m_Rotation[0]); + m_Controls->m_YrotBox->setValue(parameters.m_FiberGen.m_Rotation[1]); + m_Controls->m_ZrotBox->setValue(parameters.m_FiberGen.m_Rotation[2]); + m_Controls->m_XtransBox->setValue(parameters.m_FiberGen.m_Translation[0]); + m_Controls->m_YtransBox->setValue(parameters.m_FiberGen.m_Translation[1]); + m_Controls->m_ZtransBox->setValue(parameters.m_FiberGen.m_Translation[2]); + m_Controls->m_XscaleBox->setValue(parameters.m_FiberGen.m_Scale[0]); + m_Controls->m_YscaleBox->setValue(parameters.m_FiberGen.m_Scale[1]); + m_Controls->m_ZscaleBox->setValue(parameters.m_FiberGen.m_Scale[2]); // image generation parameters - m_Controls->m_SizeX->setValue(ffParamaters.m_ImageRegion.GetSize(0)); - m_Controls->m_SizeY->setValue(ffParamaters.m_ImageRegion.GetSize(1)); - m_Controls->m_SizeZ->setValue(ffParamaters.m_ImageRegion.GetSize(2)); - m_Controls->m_SpacingX->setValue(ffParamaters.m_ImageSpacing[0]); - m_Controls->m_SpacingY->setValue(ffParamaters.m_ImageSpacing[1]); - m_Controls->m_SpacingZ->setValue(ffParamaters.m_ImageSpacing[2]); - m_Controls->m_NumGradientsBox->setValue(ffParamaters.GetNumWeightedVolumes()); - m_Controls->m_BvalueBox->setValue(ffParamaters.m_Bvalue); - m_Controls->m_AdvancedOptionsBox_2->setChecked(ffParamaters.m_AdvancedOptions); - m_Controls->m_SignalScaleBox->setValue(ffParamaters.m_SignalScale); - m_Controls->m_TEbox->setValue(ffParamaters.m_tEcho); - m_Controls->m_LineReadoutTimeBox->setValue(ffParamaters.m_tLine); - m_Controls->m_T2starBox->setValue(ffParamaters.m_tInhom); - m_Controls->m_FiberRadius->setValue(ffParamaters.m_AxonRadius); - m_Controls->m_RelaxationBox->setChecked(ffParamaters.m_DoSimulateRelaxation); - m_Controls->m_EnforcePureFiberVoxelsBox->setChecked(ffParamaters.m_DoDisablePartialVolume); - m_Controls->m_VolumeFractionsBox->setChecked(ffParamaters.m_OutputVolumeFractions); - - if (ffParamaters.m_NoiseModel!=NULL) - { - m_Controls->m_AddNoise->setChecked(ffParamaters.m_DoAddNoise); - if (dynamic_cast*>(ffParamaters.m_NoiseModel)) + m_Controls->m_SizeX->setValue(parameters.m_SignalGen.m_ImageRegion.GetSize(0)); + m_Controls->m_SizeY->setValue(parameters.m_SignalGen.m_ImageRegion.GetSize(1)); + m_Controls->m_SizeZ->setValue(parameters.m_SignalGen.m_ImageRegion.GetSize(2)); + m_Controls->m_SpacingX->setValue(parameters.m_SignalGen.m_ImageSpacing[0]); + m_Controls->m_SpacingY->setValue(parameters.m_SignalGen.m_ImageSpacing[1]); + m_Controls->m_SpacingZ->setValue(parameters.m_SignalGen.m_ImageSpacing[2]); + m_Controls->m_NumGradientsBox->setValue(parameters.GetNumWeightedVolumes()); + m_Controls->m_BvalueBox->setValue(parameters.m_SignalGen.m_Bvalue); + m_Controls->m_AdvancedOptionsBox_2->setChecked(parameters.m_Misc.m_AdvancedOptions); + m_Controls->m_SignalScaleBox->setValue(parameters.m_SignalGen.m_SignalScale); + m_Controls->m_TEbox->setValue(parameters.m_SignalGen.m_tEcho); + m_Controls->m_LineReadoutTimeBox->setValue(parameters.m_SignalGen.m_tLine); + m_Controls->m_T2starBox->setValue(parameters.m_SignalGen.m_tInhom); + m_Controls->m_FiberRadius->setValue(parameters.m_SignalGen.m_AxonRadius); + m_Controls->m_RelaxationBox->setChecked(parameters.m_SignalGen.m_DoSimulateRelaxation); + m_Controls->m_EnforcePureFiberVoxelsBox->setChecked(parameters.m_SignalGen.m_DoDisablePartialVolume); + m_Controls->m_VolumeFractionsBox->setChecked(parameters.m_Misc.m_OutputVolumeFractions); + + if (parameters.m_NoiseModel!=NULL) + { + m_Controls->m_AddNoise->setChecked(parameters.m_SignalGen.m_DoAddNoise); + if (dynamic_cast*>(parameters.m_NoiseModel)) m_Controls->m_NoiseDistributionBox->setCurrentIndex(0); - else if (dynamic_cast*>(ffParamaters.m_NoiseModel)) + else if (dynamic_cast*>(parameters.m_NoiseModel)) m_Controls->m_NoiseDistributionBox->setCurrentIndex(1); - m_Controls->m_NoiseLevel->setValue(ffParamaters.m_NoiseModel->GetNoiseVariance()); + m_Controls->m_NoiseLevel->setValue(parameters.m_NoiseModel->GetNoiseVariance()); } else m_Controls->m_AddNoise->setChecked(false); - m_Controls->m_AddGhosts->setChecked(ffParamaters.m_DoAddGhosts); - m_Controls->m_kOffsetBox->setValue(ffParamaters.m_KspaceLineOffset); - m_Controls->m_AddAliasing->setChecked(ffParamaters.m_DoAddAliasing); - m_Controls->m_WrapBox->setValue(100*(1-ffParamaters.m_CroppingFactor)); - m_Controls->m_AddDistortions->setChecked(ffParamaters.m_DoAddDistortions); - m_Controls->m_AddSpikes->setChecked(ffParamaters.m_DoAddSpikes); - m_Controls->m_SpikeNumBox->setValue(ffParamaters.m_Spikes); - m_Controls->m_SpikeScaleBox->setValue(ffParamaters.m_SpikeAmplitude); - m_Controls->m_AddEddy->setChecked(ffParamaters.m_DoAddEddyCurrents); - m_Controls->m_EddyGradientStrength->setValue(ffParamaters.m_EddyStrength); - m_Controls->m_AddGibbsRinging->setChecked(ffParamaters.m_DoAddGibbsRinging); - m_Controls->m_AddMotion->setChecked(ffParamaters.m_DoAddMotion); - m_Controls->m_RandomMotion->setChecked(ffParamaters.m_DoRandomizeMotion); - m_Controls->m_MaxTranslationBoxX->setValue(ffParamaters.m_Translation[0]); - m_Controls->m_MaxTranslationBoxY->setValue(ffParamaters.m_Translation[1]); - m_Controls->m_MaxTranslationBoxZ->setValue(ffParamaters.m_Translation[2]); - m_Controls->m_MaxRotationBoxX->setValue(ffParamaters.m_Rotation[0]); - m_Controls->m_MaxRotationBoxY->setValue(ffParamaters.m_Rotation[1]); - m_Controls->m_MaxRotationBoxZ->setValue(ffParamaters.m_Rotation[2]); - m_Controls->m_DiffusionDirectionBox->setCurrentIndex(ffParamaters.m_DiffusionDirectionMode); - m_Controls->m_SeparationAngleBox->setValue(ffParamaters.m_FiberSeparationThreshold); + m_Controls->m_AddGhosts->setChecked(parameters.m_SignalGen.m_DoAddGhosts); + m_Controls->m_kOffsetBox->setValue(parameters.m_SignalGen.m_KspaceLineOffset); + m_Controls->m_AddAliasing->setChecked(parameters.m_SignalGen.m_DoAddAliasing); + m_Controls->m_WrapBox->setValue(100*(1-parameters.m_SignalGen.m_CroppingFactor)); + m_Controls->m_AddDistortions->setChecked(parameters.m_SignalGen.m_DoAddDistortions); + m_Controls->m_AddSpikes->setChecked(parameters.m_SignalGen.m_DoAddSpikes); + m_Controls->m_SpikeNumBox->setValue(parameters.m_SignalGen.m_Spikes); + m_Controls->m_SpikeScaleBox->setValue(parameters.m_SignalGen.m_SpikeAmplitude); + m_Controls->m_AddEddy->setChecked(parameters.m_SignalGen.m_DoAddEddyCurrents); + m_Controls->m_EddyGradientStrength->setValue(parameters.m_SignalGen.m_EddyStrength); + m_Controls->m_AddGibbsRinging->setChecked(parameters.m_SignalGen.m_DoAddGibbsRinging); + m_Controls->m_AddMotion->setChecked(parameters.m_SignalGen.m_DoAddMotion); + m_Controls->m_RandomMotion->setChecked(parameters.m_SignalGen.m_DoRandomizeMotion); + m_Controls->m_MaxTranslationBoxX->setValue(parameters.m_SignalGen.m_Translation[0]); + m_Controls->m_MaxTranslationBoxY->setValue(parameters.m_SignalGen.m_Translation[1]); + m_Controls->m_MaxTranslationBoxZ->setValue(parameters.m_SignalGen.m_Translation[2]); + m_Controls->m_MaxRotationBoxX->setValue(parameters.m_SignalGen.m_Rotation[0]); + m_Controls->m_MaxRotationBoxY->setValue(parameters.m_SignalGen.m_Rotation[1]); + m_Controls->m_MaxRotationBoxZ->setValue(parameters.m_SignalGen.m_Rotation[2]); + m_Controls->m_DiffusionDirectionBox->setCurrentIndex(parameters.m_SignalGen.m_DiffusionDirectionMode); + m_Controls->m_SeparationAngleBox->setValue(parameters.m_SignalGen.m_FiberSeparationThreshold); m_Controls->m_Compartment1Box->setCurrentIndex(0); m_Controls->m_Compartment2Box->setCurrentIndex(0); m_Controls->m_Compartment3Box->setCurrentIndex(0); m_Controls->m_Compartment4Box->setCurrentIndex(0); - for (unsigned int i=0; i* signalModel = NULL; - if (im_CompartmentId) { case 1: { if (dynamic_cast*>(signalModel)) { mitk::StickModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_StickWidget1->SetT2(model->GetT2()); m_Controls->m_StickWidget1->SetD(model->GetDiffusivity()); m_Controls->m_Compartment1Box->setCurrentIndex(0); break; } else if (dynamic_cast*>(signalModel)) { mitk::TensorModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_TensorWidget1->SetT2(model->GetT2()); m_Controls->m_TensorWidget1->SetD1(model->GetDiffusivity1()); m_Controls->m_TensorWidget1->SetD2(model->GetDiffusivity1()); m_Controls->m_TensorWidget1->SetD3(model->GetDiffusivity1()); m_Controls->m_Compartment1Box->setCurrentIndex(2); break; } else if (dynamic_cast*>(signalModel)) { mitk::RawShModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_PrototypeWidget1->SetNumberOfSamples(model->GetMaxNumKernels()); m_Controls->m_PrototypeWidget1->SetMinFa(model->GetFaRange().first); m_Controls->m_PrototypeWidget1->SetMaxFa(model->GetFaRange().second); m_Controls->m_PrototypeWidget1->SetMinAdc(model->GetAdcRange().first); m_Controls->m_PrototypeWidget1->SetMaxAdc(model->GetAdcRange().second); m_Controls->m_Compartment1Box->setCurrentIndex(3); break; } break; } case 2: { if (dynamic_cast*>(signalModel)) { mitk::StickModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_StickWidget2->SetT2(model->GetT2()); m_Controls->m_StickWidget2->SetD(model->GetDiffusivity()); m_Controls->m_Compartment2Box->setCurrentIndex(1); break; } else if (dynamic_cast*>(signalModel)) { mitk::TensorModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_TensorWidget2->SetT2(model->GetT2()); m_Controls->m_TensorWidget2->SetD1(model->GetDiffusivity1()); m_Controls->m_TensorWidget2->SetD2(model->GetDiffusivity1()); m_Controls->m_TensorWidget2->SetD3(model->GetDiffusivity1()); m_Controls->m_Compartment2Box->setCurrentIndex(3); break; } break; } case 3: { if (dynamic_cast*>(signalModel)) { mitk::BallModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_BallWidget1->SetT2(model->GetT2()); m_Controls->m_BallWidget1->SetD(model->GetDiffusivity()); m_Controls->m_Compartment3Box->setCurrentIndex(0); break; } else if (dynamic_cast*>(signalModel)) { mitk::AstroStickModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_AstrosticksWidget1->SetT2(model->GetT2()); m_Controls->m_AstrosticksWidget1->SetD(model->GetDiffusivity()); m_Controls->m_AstrosticksWidget1->SetRandomizeSticks(model->GetRandomizeSticks()); m_Controls->m_Compartment3Box->setCurrentIndex(1); break; } else if (dynamic_cast*>(signalModel)) { mitk::DotModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_DotWidget1->SetT2(model->GetT2()); m_Controls->m_Compartment3Box->setCurrentIndex(2); break; } else if (dynamic_cast*>(signalModel)) { mitk::RawShModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_PrototypeWidget3->SetNumberOfSamples(model->GetMaxNumKernels()); m_Controls->m_PrototypeWidget3->SetMinFa(model->GetFaRange().first); m_Controls->m_PrototypeWidget3->SetMaxFa(model->GetFaRange().second); m_Controls->m_PrototypeWidget3->SetMinAdc(model->GetAdcRange().first); m_Controls->m_PrototypeWidget3->SetMaxAdc(model->GetAdcRange().second); m_Controls->m_Compartment3Box->setCurrentIndex(3); break; } break; } case 4: { if (dynamic_cast*>(signalModel)) { mitk::BallModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_BallWidget2->SetT2(model->GetT2()); m_Controls->m_BallWidget2->SetD(model->GetDiffusivity()); m_Controls->m_Compartment4Box->setCurrentIndex(1); break; } else if (dynamic_cast*>(signalModel)) { mitk::AstroStickModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_AstrosticksWidget2->SetT2(model->GetT2()); m_Controls->m_AstrosticksWidget2->SetD(model->GetDiffusivity()); m_Controls->m_AstrosticksWidget2->SetRandomizeSticks(model->GetRandomizeSticks()); m_Controls->m_Compartment4Box->setCurrentIndex(2); break; } else if (dynamic_cast*>(signalModel)) { mitk::DotModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_DotWidget2->SetT2(model->GetT2()); m_Controls->m_Compartment4Box->setCurrentIndex(3); break; } else if (dynamic_cast*>(signalModel)) { mitk::RawShModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_PrototypeWidget4->SetNumberOfSamples(model->GetMaxNumKernels()); m_Controls->m_PrototypeWidget4->SetMinFa(model->GetFaRange().first); m_Controls->m_PrototypeWidget4->SetMaxFa(model->GetFaRange().second); m_Controls->m_PrototypeWidget4->SetMinAdc(model->GetAdcRange().first); m_Controls->m_PrototypeWidget4->SetMaxAdc(model->GetAdcRange().second); m_Controls->m_Compartment4Box->setCurrentIndex(4); break; } break; } } } } void QmitkFiberfoxView::ShowAdvancedOptions(int state) { if (state) { m_Controls->m_AdvancedFiberOptionsFrame->setVisible(true); m_Controls->m_AdvancedSignalOptionsFrame->setVisible(true); m_Controls->m_AdvancedOptionsBox->setChecked(true); m_Controls->m_AdvancedOptionsBox_2->setChecked(true); } else { m_Controls->m_AdvancedFiberOptionsFrame->setVisible(false); m_Controls->m_AdvancedSignalOptionsFrame->setVisible(false); m_Controls->m_AdvancedOptionsBox->setChecked(false); m_Controls->m_AdvancedOptionsBox_2->setChecked(false); } } void QmitkFiberfoxView::Comp1ModelFrameVisibility(int index) { m_Controls->m_StickWidget1->setVisible(false); m_Controls->m_ZeppelinWidget1->setVisible(false); m_Controls->m_TensorWidget1->setVisible(false); m_Controls->m_PrototypeWidget1->setVisible(false); switch (index) { case 0: m_Controls->m_StickWidget1->setVisible(true); break; case 1: m_Controls->m_ZeppelinWidget1->setVisible(true); break; case 2: m_Controls->m_TensorWidget1->setVisible(true); break; case 3: m_Controls->m_PrototypeWidget1->setVisible(true); break; } } void QmitkFiberfoxView::Comp2ModelFrameVisibility(int index) { m_Controls->m_StickWidget2->setVisible(false); m_Controls->m_ZeppelinWidget2->setVisible(false); m_Controls->m_TensorWidget2->setVisible(false); switch (index) { case 0: break; case 1: m_Controls->m_StickWidget2->setVisible(true); break; case 2: m_Controls->m_ZeppelinWidget2->setVisible(true); break; case 3: m_Controls->m_TensorWidget2->setVisible(true); break; } } void QmitkFiberfoxView::Comp3ModelFrameVisibility(int index) { m_Controls->m_BallWidget1->setVisible(false); m_Controls->m_AstrosticksWidget1->setVisible(false); m_Controls->m_DotWidget1->setVisible(false); m_Controls->m_PrototypeWidget3->setVisible(false); switch (index) { case 0: m_Controls->m_BallWidget1->setVisible(true); break; case 1: m_Controls->m_AstrosticksWidget1->setVisible(true); break; case 2: m_Controls->m_DotWidget1->setVisible(true); break; case 3: m_Controls->m_PrototypeWidget3->setVisible(true); break; } } void QmitkFiberfoxView::Comp4ModelFrameVisibility(int index) { m_Controls->m_BallWidget2->setVisible(false); m_Controls->m_AstrosticksWidget2->setVisible(false); m_Controls->m_DotWidget2->setVisible(false); m_Controls->m_PrototypeWidget4->setVisible(false); m_Controls->m_Comp4FractionFrame->setVisible(false); switch (index) { case 0: break; case 1: m_Controls->m_BallWidget2->setVisible(true); m_Controls->m_Comp4FractionFrame->setVisible(true); break; case 2: m_Controls->m_AstrosticksWidget2->setVisible(true); m_Controls->m_Comp4FractionFrame->setVisible(true); break; case 3: m_Controls->m_DotWidget2->setVisible(true); m_Controls->m_Comp4FractionFrame->setVisible(true); break; case 4: m_Controls->m_PrototypeWidget4->setVisible(true); m_Controls->m_Comp4FractionFrame->setVisible(true); break; } } void QmitkFiberfoxView::OnConstantRadius(int value) { if (value>0 && m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnAddMotion(int value) { if (value>0) m_Controls->m_MotionArtifactFrame->setVisible(true); else m_Controls->m_MotionArtifactFrame->setVisible(false); } void QmitkFiberfoxView::OnAddAliasing(int value) { if (value>0) m_Controls->m_AliasingFrame->setVisible(true); else m_Controls->m_AliasingFrame->setVisible(false); } void QmitkFiberfoxView::OnAddSpikes(int value) { if (value>0) m_Controls->m_SpikeFrame->setVisible(true); else m_Controls->m_SpikeFrame->setVisible(false); } void QmitkFiberfoxView::OnAddEddy(int value) { if (value>0) m_Controls->m_EddyFrame->setVisible(true); else m_Controls->m_EddyFrame->setVisible(false); } void QmitkFiberfoxView::OnAddDistortions(int value) { if (value>0) m_Controls->m_DistortionsFrame->setVisible(true); else m_Controls->m_DistortionsFrame->setVisible(false); } void QmitkFiberfoxView::OnAddGhosts(int value) { if (value>0) m_Controls->m_GhostFrame->setVisible(true); else m_Controls->m_GhostFrame->setVisible(false); } void QmitkFiberfoxView::OnAddNoise(int value) { if (value>0) m_Controls->m_NoiseFrame->setVisible(true); else m_Controls->m_NoiseFrame->setVisible(false); } void QmitkFiberfoxView::OnDistributionChanged(int value) { if (value==1) m_Controls->m_VarianceBox->setVisible(true); else m_Controls->m_VarianceBox->setVisible(false); if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnVarianceChanged(double) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnFiberDensityChanged(int) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnFiberSamplingChanged(double) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnTensionChanged(double) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnContinuityChanged(double) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnBiasChanged(double) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::AlignOnGrid() { for (unsigned int i=0; i(m_SelectedFiducials.at(i)->GetData()); mitk::Point3D wc0 = pe->GetWorldControlPoint(0); mitk::DataStorage::SetOfObjects::ConstPointer parentFibs = GetDataStorage()->GetSources(m_SelectedFiducials.at(i)); for( mitk::DataStorage::SetOfObjects::const_iterator it = parentFibs->begin(); it != parentFibs->end(); ++it ) { mitk::DataNode::Pointer pFibNode = *it; if ( pFibNode.IsNotNull() && dynamic_cast(pFibNode->GetData()) ) { mitk::DataStorage::SetOfObjects::ConstPointer parentImgs = GetDataStorage()->GetSources(pFibNode); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = parentImgs->begin(); it2 != parentImgs->end(); ++it2 ) { mitk::DataNode::Pointer pImgNode = *it2; if ( pImgNode.IsNotNull() && dynamic_cast(pImgNode->GetData()) ) { mitk::Image::Pointer img = dynamic_cast(pImgNode->GetData()); mitk::BaseGeometry::Pointer geom = img->GetGeometry(); itk::Index<3> idx; geom->WorldToIndex(wc0, idx); mitk::Point3D cIdx; cIdx[0]=idx[0]; cIdx[1]=idx[1]; cIdx[2]=idx[2]; mitk::Point3D world; geom->IndexToWorld(cIdx,world); mitk::Vector3D trans = world - wc0; pe->GetGeometry()->Translate(trans); break; } } break; } } } for(unsigned int i=0; iGetSources(fibNode); for( mitk::DataStorage::SetOfObjects::const_iterator it = sources->begin(); it != sources->end(); ++it ) { mitk::DataNode::Pointer imgNode = *it; if ( imgNode.IsNotNull() && dynamic_cast(imgNode->GetData()) ) { mitk::DataStorage::SetOfObjects::ConstPointer derivations = GetDataStorage()->GetDerivations(fibNode); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = derivations->begin(); it2 != derivations->end(); ++it2 ) { mitk::DataNode::Pointer fiducialNode = *it2; if ( fiducialNode.IsNotNull() && dynamic_cast(fiducialNode->GetData()) ) { mitk::PlanarEllipse::Pointer pe = dynamic_cast(fiducialNode->GetData()); mitk::Point3D wc0 = pe->GetWorldControlPoint(0); mitk::Image::Pointer img = dynamic_cast(imgNode->GetData()); mitk::BaseGeometry::Pointer geom = img->GetGeometry(); itk::Index<3> idx; geom->WorldToIndex(wc0, idx); mitk::Point3D cIdx; cIdx[0]=idx[0]; cIdx[1]=idx[1]; cIdx[2]=idx[2]; mitk::Point3D world; geom->IndexToWorld(cIdx,world); mitk::Vector3D trans = world - wc0; pe->GetGeometry()->Translate(trans); } } break; } } } for(unsigned int i=0; i(m_SelectedImages.at(i)->GetData()); mitk::DataStorage::SetOfObjects::ConstPointer derivations = GetDataStorage()->GetDerivations(m_SelectedImages.at(i)); for( mitk::DataStorage::SetOfObjects::const_iterator it = derivations->begin(); it != derivations->end(); ++it ) { mitk::DataNode::Pointer fibNode = *it; if ( fibNode.IsNotNull() && dynamic_cast(fibNode->GetData()) ) { mitk::DataStorage::SetOfObjects::ConstPointer derivations2 = GetDataStorage()->GetDerivations(fibNode); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = derivations2->begin(); it2 != derivations2->end(); ++it2 ) { mitk::DataNode::Pointer fiducialNode = *it2; if ( fiducialNode.IsNotNull() && dynamic_cast(fiducialNode->GetData()) ) { mitk::PlanarEllipse::Pointer pe = dynamic_cast(fiducialNode->GetData()); mitk::Point3D wc0 = pe->GetWorldControlPoint(0); mitk::BaseGeometry::Pointer geom = img->GetGeometry(); itk::Index<3> idx; geom->WorldToIndex(wc0, idx); mitk::Point3D cIdx; cIdx[0]=idx[0]; cIdx[1]=idx[1]; cIdx[2]=idx[2]; mitk::Point3D world; geom->IndexToWorld(cIdx,world); mitk::Vector3D trans = world - wc0; pe->GetGeometry()->Translate(trans); } } } } } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnFlipButton() { if (m_SelectedFiducial.IsNull()) return; std::map::iterator it = m_DataNodeToPlanarFigureData.find(m_SelectedFiducial.GetPointer()); if( it != m_DataNodeToPlanarFigureData.end() ) { QmitkPlanarFigureData& data = it->second; data.m_Flipped += 1; data.m_Flipped %= 2; } if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } QmitkFiberfoxView::GradientListType QmitkFiberfoxView::GenerateHalfShell(int NPoints) { NPoints *= 2; GradientListType pointshell; int numB0 = NPoints/20; if (numB0==0) numB0=1; GradientType g; g.Fill(0.0); for (int i=0; i theta; theta.set_size(NPoints); vnl_vector phi; phi.set_size(NPoints); double C = sqrt(4*M_PI); phi(0) = 0.0; phi(NPoints-1) = 0.0; for(int i=0; i0 && i std::vector > QmitkFiberfoxView::MakeGradientList() { std::vector > retval; vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); // Add 0 vector for B0 int numB0 = ndirs/10; if (numB0==0) numB0=1; itk::Vector v; v.Fill(0.0); for (int i=0; i v; v[0] = U->get(0,i); v[1] = U->get(1,i); v[2] = U->get(2,i); retval.push_back(v); } return retval; } void QmitkFiberfoxView::OnAddBundle() { if (m_SelectedImage.IsNull()) return; mitk::DataStorage::SetOfObjects::ConstPointer children = GetDataStorage()->GetDerivations(m_SelectedImage); mitk::FiberBundleX::Pointer bundle = mitk::FiberBundleX::New(); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( bundle ); QString name = QString("Bundle_%1").arg(children->size()); node->SetName(name.toStdString()); m_SelectedBundles.push_back(node); UpdateGui(); GetDataStorage()->Add(node, m_SelectedImage); } void QmitkFiberfoxView::OnDrawROI() { if (m_SelectedBundles.empty()) OnAddBundle(); if (m_SelectedBundles.empty()) return; mitk::DataStorage::SetOfObjects::ConstPointer children = GetDataStorage()->GetDerivations(m_SelectedBundles.at(0)); mitk::PlanarEllipse::Pointer figure = mitk::PlanarEllipse::New(); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( figure ); node->SetBoolProperty("planarfigure.3drendering", true); QList nodes = this->GetDataManagerSelection(); for( int i=0; iSetSelected(false); m_SelectedFiducial = node; QString name = QString("Fiducial_%1").arg(children->size()); node->SetName(name.toStdString()); node->SetSelected(true); this->DisableCrosshairNavigation(); mitk::PlanarFigureInteractor::Pointer figureInteractor = dynamic_cast(node->GetDataInteractor().GetPointer()); if(figureInteractor.IsNull()) { figureInteractor = mitk::PlanarFigureInteractor::New(); us::Module* planarFigureModule = us::ModuleRegistry::GetModule( "MitkPlanarFigure" ); figureInteractor->LoadStateMachine("PlanarFigureInteraction.xml", planarFigureModule ); figureInteractor->SetEventConfig( "PlanarFigureConfig.xml", planarFigureModule ); figureInteractor->SetDataNode( node ); } UpdateGui(); GetDataStorage()->Add(node, m_SelectedBundles.at(0)); } bool CompareLayer(mitk::DataNode::Pointer i,mitk::DataNode::Pointer j) { int li = -1; i->GetPropertyValue("layer", li); int lj = -1; j->GetPropertyValue("layer", lj); return liGetSources(m_SelectedFiducial); for( mitk::DataStorage::SetOfObjects::const_iterator it = parents->begin(); it != parents->end(); ++it ) if(dynamic_cast((*it)->GetData())) m_SelectedBundles.push_back(*it); if (m_SelectedBundles.empty()) return; } vector< vector< mitk::PlanarEllipse::Pointer > > fiducials; vector< vector< unsigned int > > fliplist; for (unsigned int i=0; iGetDerivations(m_SelectedBundles.at(i)); std::vector< mitk::DataNode::Pointer > childVector; for( mitk::DataStorage::SetOfObjects::const_iterator it = children->begin(); it != children->end(); ++it ) childVector.push_back(*it); sort(childVector.begin(), childVector.end(), CompareLayer); vector< mitk::PlanarEllipse::Pointer > fib; vector< unsigned int > flip; float radius = 1; int count = 0; for( std::vector< mitk::DataNode::Pointer >::const_iterator it = childVector.begin(); it != childVector.end(); ++it ) { mitk::DataNode::Pointer node = *it; if ( node.IsNotNull() && dynamic_cast(node->GetData()) ) { mitk::PlanarEllipse* ellipse = dynamic_cast(node->GetData()); if (m_Controls->m_ConstantRadiusBox->isChecked()) { ellipse->SetTreatAsCircle(true); mitk::Point2D c = ellipse->GetControlPoint(0); mitk::Point2D p = ellipse->GetControlPoint(1); mitk::Vector2D v = p-c; if (count==0) { radius = v.GetVnlVector().magnitude(); ellipse->SetControlPoint(1, p); } else { v.Normalize(); v *= radius; ellipse->SetControlPoint(1, c+v); } } fib.push_back(ellipse); std::map::iterator it = m_DataNodeToPlanarFigureData.find(node.GetPointer()); if( it != m_DataNodeToPlanarFigureData.end() ) { QmitkPlanarFigureData& data = it->second; flip.push_back(data.m_Flipped); } else flip.push_back(0); } count++; } if (fib.size()>1) { fiducials.push_back(fib); fliplist.push_back(flip); } else if (fib.size()>0) m_SelectedBundles.at(i)->SetData( mitk::FiberBundleX::New() ); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } itk::FibersFromPlanarFiguresFilter::Pointer filter = itk::FibersFromPlanarFiguresFilter::New(); filter->SetFiducials(fiducials); filter->SetFlipList(fliplist); switch(m_Controls->m_DistributionBox->currentIndex()){ case 0: filter->SetFiberDistribution(itk::FibersFromPlanarFiguresFilter::DISTRIBUTE_UNIFORM); break; case 1: filter->SetFiberDistribution(itk::FibersFromPlanarFiguresFilter::DISTRIBUTE_GAUSSIAN); filter->SetVariance(m_Controls->m_VarianceBox->value()); break; } filter->SetDensity(m_Controls->m_FiberDensityBox->value()); filter->SetTension(m_Controls->m_TensionBox->value()); filter->SetContinuity(m_Controls->m_ContinuityBox->value()); filter->SetBias(m_Controls->m_BiasBox->value()); filter->SetFiberSampling(m_Controls->m_FiberSamplingBox->value()); filter->Update(); vector< mitk::FiberBundleX::Pointer > fiberBundles = filter->GetFiberBundles(); for (unsigned int i=0; iSetData( fiberBundles.at(i) ); if (fiberBundles.at(i)->GetNumFibers()>50000) m_SelectedBundles.at(i)->SetVisibility(false); } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberfoxView::GenerateImage() { if (m_SelectedBundles.empty() && m_SelectedDWI.IsNull()) { mitk::Image::Pointer image = mitk::ImageGenerator::GenerateGradientImage( m_Controls->m_SizeX->value(), m_Controls->m_SizeY->value(), m_Controls->m_SizeZ->value(), m_Controls->m_SpacingX->value(), m_Controls->m_SpacingY->value(), m_Controls->m_SpacingZ->value()); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); node->SetName("Dummy"); unsigned int window = m_Controls->m_SizeX->value()*m_Controls->m_SizeY->value()*m_Controls->m_SizeZ->value(); unsigned int level = window/2; mitk::LevelWindow lw; lw.SetLevelWindow(level, window); node->SetProperty( "levelwindow", mitk::LevelWindowProperty::New( lw ) ); GetDataStorage()->Add(node); m_SelectedImage = 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(); } else if (!m_SelectedBundles.empty()) SimulateImageFromFibers(m_SelectedBundles.at(0)); else if (m_SelectedDWI.IsNotNull()) SimulateForExistingDwi(m_SelectedDWI); } void QmitkFiberfoxView::SimulateForExistingDwi(mitk::DataNode* imageNode) { if (!dynamic_cast*>(imageNode->GetData())) return; FiberfoxParameters parameters = UpdateImageParameters(); if (parameters.m_NoiseModel==NULL && - parameters.m_Spikes==0 && - parameters.m_FrequencyMap.IsNull() && - parameters.m_KspaceLineOffset<=0.000001 && - !parameters.m_DoAddGibbsRinging && - !(parameters.m_EddyStrength>0) && - parameters.m_CroppingFactor>0.999) + parameters.m_SignalGen.m_Spikes==0 && + parameters.m_SignalGen.m_FrequencyMap.IsNull() && + parameters.m_SignalGen.m_KspaceLineOffset<=0.000001 && + !parameters.m_SignalGen.m_DoAddGibbsRinging && + !(parameters.m_SignalGen.m_EddyStrength>0) && + parameters.m_SignalGen.m_CroppingFactor>0.999) { QMessageBox::information( NULL, "Simulation cancelled", "No valid artifact enabled! Motion artifacts and relaxation effects can NOT be added to an existing diffusion weighted image."); return; } mitk::DiffusionImage::Pointer diffImg = dynamic_cast*>(imageNode->GetData()); m_ArtifactsToDwiFilter = itk::AddArtifactsToDwiImageFilter< short >::New(); m_ArtifactsToDwiFilter->SetInput(diffImg->GetVectorImage()); - parameters.m_ParentNode = imageNode; + parameters.m_Misc.m_ParentNode = imageNode; m_ArtifactsToDwiFilter->SetParameters(parameters); m_Worker.m_FilterType = 1; m_Thread.start(QThread::LowestPriority); } void QmitkFiberfoxView::SimulateImageFromFibers(mitk::DataNode* fiberNode) { mitk::FiberBundleX::Pointer fiberBundle = dynamic_cast(fiberNode->GetData()); if (fiberBundle->GetNumFibers()<=0) return; FiberfoxParameters parameters = UpdateImageParameters(); m_TractsToDwiFilter = itk::TractsToDWIImageFilter< short >::New(); - parameters.m_ParentNode = fiberNode; + parameters.m_Misc.m_ParentNode = fiberNode; if (m_SelectedDWI.IsNotNull()) { mitk::DiffusionImage::Pointer diffImg = dynamic_cast*>(m_SelectedDWI->GetData()); bool doSampling = false; for (unsigned int i=0; i* >(parameters.m_FiberModelList[i]) ) doSampling = true; for (unsigned int i=0; i* >(parameters.m_NonFiberModelList[i]) ) doSampling = true; // sample prototype signals if ( doSampling ) { const int shOrder = 2; typedef itk::DiffusionTensor3DReconstructionImageFilter< short, short, double > TensorReconstructionImageFilterType; TensorReconstructionImageFilterType::Pointer filter = TensorReconstructionImageFilterType::New(); filter->SetGradientImage( diffImg->GetDirections(), diffImg->GetVectorImage() ); filter->SetBValue(diffImg->GetReferenceBValue()); filter->Update(); itk::Image< itk::DiffusionTensor3D< double >, 3 >::Pointer tensorImage = filter->GetOutput(); const int NumCoeffs = (shOrder*shOrder + shOrder + 2)/2 + shOrder; typedef itk::AnalyticalDiffusionQballReconstructionImageFilter QballFilterType; QballFilterType::Pointer qballfilter = QballFilterType::New(); qballfilter->SetGradientImage( diffImg->GetDirections(), diffImg->GetVectorImage() ); qballfilter->SetBValue(diffImg->GetReferenceBValue()); qballfilter->SetLambda(0.006); qballfilter->SetNormalizationMethod(QballFilterType::QBAR_RAW_SIGNAL); qballfilter->Update(); QballFilterType::CoefficientImageType::Pointer itkFeatureImage = qballfilter->GetCoefficientImage(); itk::AdcImageFilter< short, double >::Pointer adcFilter = itk::AdcImageFilter< short, double >::New(); adcFilter->SetInput(diffImg->GetVectorImage()); adcFilter->SetGradientDirections(diffImg->GetDirections()); adcFilter->SetB_value(diffImg->GetReferenceBValue()); adcFilter->Update(); ItkDoubleImgType::Pointer adcImage = adcFilter->GetOutput(); int b0Index; for (unsigned int i=0; iGetDirectionsWithoutMeasurementFrame()->Size(); i++) if ( diffImg->GetDirectionsWithoutMeasurementFrame()->GetElement(i).magnitude()<0.001 ) { b0Index = i; break; } double max = 0; { itk::ImageRegionIterator< itk::VectorImage< short, 3 > > it(diffImg->GetVectorImage(), diffImg->GetVectorImage()->GetLargestPossibleRegion()); while(!it.IsAtEnd()) { - if (parameters.m_MaskImage.IsNotNull() && parameters.m_MaskImage->GetPixel(it.GetIndex())<=0) + if (parameters.m_SignalGen.m_MaskImage.IsNotNull() && parameters.m_SignalGen.m_MaskImage->GetPixel(it.GetIndex())<=0) { ++it; continue; } if (it.Get()[b0Index]>max) max = it.Get()[b0Index]; ++it; } } MITK_INFO << "Sampling signal kernels."; itk::ImageRegionIterator< itk::Image< itk::DiffusionTensor3D< double >, 3 > > it(tensorImage, tensorImage->GetLargestPossibleRegion()); while(!it.IsAtEnd()) { bool skipPixel = false; - if (parameters.m_MaskImage.IsNotNull() && parameters.m_MaskImage->GetPixel(it.GetIndex())<=0) + if (parameters.m_SignalGen.m_MaskImage.IsNotNull() && parameters.m_SignalGen.m_MaskImage->GetPixel(it.GetIndex())<=0) { ++it; continue; } for (unsigned int i=0; iGetDirections()->Size(); i++) { if (diffImg->GetVectorImage()->GetPixel(it.GetIndex())[i]!=diffImg->GetVectorImage()->GetPixel(it.GetIndex())[i] || diffImg->GetVectorImage()->GetPixel(it.GetIndex())[i]<=0 || diffImg->GetVectorImage()->GetPixel(it.GetIndex())[i]>diffImg->GetVectorImage()->GetPixel(it.GetIndex())[b0Index]) { skipPixel = true; break; } } if (skipPixel) { ++it; continue; } typedef itk::DiffusionTensor3D TensorType; TensorType::EigenValuesArrayType eigenvalues; TensorType::EigenVectorsMatrixType eigenvectors; TensorType tensor = it.Get(); double FA = tensor.GetFractionalAnisotropy(); double ADC = adcImage->GetPixel(it.GetIndex()); QballFilterType::CoefficientImageType::PixelType itkv = itkFeatureImage->GetPixel(it.GetIndex()); vnl_vector_fixed< double, NumCoeffs > coeffs; for (unsigned int c=0; c* model = dynamic_cast< RawShModel* >(parameters.m_FiberModelList[i]); if ( model && model->GetMaxNumKernels()>model->GetNumberOfKernels() && FA>model->GetFaRange().first && FAGetFaRange().second && ADC>model->GetAdcRange().first && ADCGetAdcRange().second) { if (model->SetShCoefficients( coeffs, (double)diffImg->GetVectorImage()->GetPixel(it.GetIndex())[b0Index]/max )) { tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); itk::Vector dir; dir[0] = eigenvectors(2, 0); dir[1] = eigenvectors(2, 1); dir[2] = eigenvectors(2, 2); model->m_PrototypeMaxDirection.push_back(dir); MITK_INFO << "WM KERNEL: " << it.GetIndex() << " (" << model->GetNumberOfKernels() << ")"; } } } for (unsigned int i=0; i* model = dynamic_cast< RawShModel* >(parameters.m_NonFiberModelList[i]); if ( model && model->GetMaxNumKernels()>model->GetNumberOfKernels() && FA>model->GetFaRange().first && FAGetFaRange().second && ADC>model->GetAdcRange().first && ADCGetAdcRange().second) { if (model->SetShCoefficients( coeffs, (double)diffImg->GetVectorImage()->GetPixel(it.GetIndex())[b0Index]/max )) MITK_INFO << "CSF/GM KERNEL: " << it.GetIndex() << " (" << model->GetNumberOfKernels() << ")"; } } ++it; } for (unsigned int i=0; i* model = dynamic_cast< RawShModel* >(parameters.m_FiberModelList[i]); if ( model && model->GetNumberOfKernels()<=0 ) { QMessageBox::information( NULL, "Simulation cancelled", "No suitable voxels found for fiber compartment "+QString::number(i)); return; } } for (unsigned int i=0; i* model = dynamic_cast< RawShModel* >(parameters.m_NonFiberModelList[i]); if ( model && model->GetNumberOfKernels()<=0 ) { QMessageBox::information( NULL, "Simulation cancelled", "No suitable voxels found for non-fiber compartment "+QString::number(i)); return; } } } } else if ( m_Controls->m_Compartment1Box->currentIndex()==3 || m_Controls->m_Compartment3Box->currentIndex()==3 || m_Controls->m_Compartment4Box->currentIndex()==4 ) { QMessageBox::information( NULL, "Simulation cancelled", "Prototype signal but no diffusion-weighted image selected to sample signal from."); return; } m_TractsToDwiFilter->SetParameters(parameters); m_TractsToDwiFilter->SetFiberBundle(fiberBundle); m_Worker.m_FilterType = 0; m_Thread.start(QThread::LowestPriority); } void QmitkFiberfoxView::ApplyTransform() { vector< mitk::DataNode::Pointer > selectedBundles; for(unsigned int i=0; iGetDerivations(m_SelectedImages.at(i)); for( mitk::DataStorage::SetOfObjects::const_iterator it = derivations->begin(); it != derivations->end(); ++it ) { mitk::DataNode::Pointer fibNode = *it; if ( fibNode.IsNotNull() && dynamic_cast(fibNode->GetData()) ) selectedBundles.push_back(fibNode); } } if (selectedBundles.empty()) selectedBundles = m_SelectedBundles2; if (!selectedBundles.empty()) { for (std::vector::const_iterator it = selectedBundles.begin(); it!=selectedBundles.end(); ++it) { mitk::FiberBundleX::Pointer fib = dynamic_cast((*it)->GetData()); fib->RotateAroundAxis(m_Controls->m_XrotBox->value(), m_Controls->m_YrotBox->value(), m_Controls->m_ZrotBox->value()); fib->TranslateFibers(m_Controls->m_XtransBox->value(), m_Controls->m_YtransBox->value(), m_Controls->m_ZtransBox->value()); fib->ScaleFibers(m_Controls->m_XscaleBox->value(), m_Controls->m_YscaleBox->value(), m_Controls->m_ZscaleBox->value()); // handle child fiducials if (m_Controls->m_IncludeFiducials->isChecked()) { mitk::DataStorage::SetOfObjects::ConstPointer derivations = GetDataStorage()->GetDerivations(*it); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = derivations->begin(); it2 != derivations->end(); ++it2 ) { mitk::DataNode::Pointer fiducialNode = *it2; if ( fiducialNode.IsNotNull() && dynamic_cast(fiducialNode->GetData()) ) { mitk::PlanarEllipse* pe = dynamic_cast(fiducialNode->GetData()); mitk::BaseGeometry* geom = pe->GetGeometry(); // translate mitk::Vector3D world; world[0] = m_Controls->m_XtransBox->value(); world[1] = m_Controls->m_YtransBox->value(); world[2] = m_Controls->m_ZtransBox->value(); geom->Translate(world); // calculate rotation matrix double x = m_Controls->m_XrotBox->value()*M_PI/180; double y = m_Controls->m_YrotBox->value()*M_PI/180; double z = m_Controls->m_ZrotBox->value()*M_PI/180; itk::Matrix< double, 3, 3 > rotX; rotX.SetIdentity(); rotX[1][1] = cos(x); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(x); rotX[2][1] = -rotX[1][2]; itk::Matrix< double, 3, 3 > rotY; rotY.SetIdentity(); rotY[0][0] = cos(y); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(y); rotY[2][0] = -rotY[0][2]; itk::Matrix< double, 3, 3 > rotZ; rotZ.SetIdentity(); rotZ[0][0] = cos(z); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(z); rotZ[1][0] = -rotZ[0][1]; itk::Matrix< double, 3, 3 > rot = rotZ*rotY*rotX; // transform control point coordinate into geometry translation geom->SetOrigin(pe->GetWorldControlPoint(0)); mitk::Point2D cp; cp.Fill(0.0); pe->SetControlPoint(0, cp); // rotate fiducial geom->GetIndexToWorldTransform()->SetMatrix(rot*geom->GetIndexToWorldTransform()->GetMatrix()); // implicit translation mitk::Vector3D trans; trans[0] = geom->GetOrigin()[0]-fib->GetGeometry()->GetCenter()[0]; trans[1] = geom->GetOrigin()[1]-fib->GetGeometry()->GetCenter()[1]; trans[2] = geom->GetOrigin()[2]-fib->GetGeometry()->GetCenter()[2]; mitk::Vector3D newWc = rot*trans; newWc = newWc-trans; geom->Translate(newWc); pe->Modified(); } } } } } else { for (unsigned int i=0; i(m_SelectedFiducials.at(i)->GetData()); mitk::BaseGeometry* geom = pe->GetGeometry(); // translate mitk::Vector3D world; world[0] = m_Controls->m_XtransBox->value(); world[1] = m_Controls->m_YtransBox->value(); world[2] = m_Controls->m_ZtransBox->value(); geom->Translate(world); // calculate rotation matrix double x = m_Controls->m_XrotBox->value()*M_PI/180; double y = m_Controls->m_YrotBox->value()*M_PI/180; double z = m_Controls->m_ZrotBox->value()*M_PI/180; itk::Matrix< double, 3, 3 > rotX; rotX.SetIdentity(); rotX[1][1] = cos(x); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(x); rotX[2][1] = -rotX[1][2]; itk::Matrix< double, 3, 3 > rotY; rotY.SetIdentity(); rotY[0][0] = cos(y); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(y); rotY[2][0] = -rotY[0][2]; itk::Matrix< double, 3, 3 > rotZ; rotZ.SetIdentity(); rotZ[0][0] = cos(z); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(z); rotZ[1][0] = -rotZ[0][1]; itk::Matrix< double, 3, 3 > rot = rotZ*rotY*rotX; // transform control point coordinate into geometry translation geom->SetOrigin(pe->GetWorldControlPoint(0)); mitk::Point2D cp; cp.Fill(0.0); pe->SetControlPoint(0, cp); // rotate fiducial geom->GetIndexToWorldTransform()->SetMatrix(rot*geom->GetIndexToWorldTransform()->GetMatrix()); pe->Modified(); } if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberfoxView::CopyBundles() { if ( m_SelectedBundles.size()<1 ){ QMessageBox::information( NULL, "Warning", "Select at least one fiber bundle!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least one fiber bundle!"; return; } for (std::vector::const_iterator it = m_SelectedBundles.begin(); it!=m_SelectedBundles.end(); ++it) { // find parent image mitk::DataNode::Pointer parentNode; mitk::DataStorage::SetOfObjects::ConstPointer parentImgs = GetDataStorage()->GetSources(*it); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = parentImgs->begin(); it2 != parentImgs->end(); ++it2 ) { mitk::DataNode::Pointer pImgNode = *it2; if ( pImgNode.IsNotNull() && dynamic_cast(pImgNode->GetData()) ) { parentNode = pImgNode; break; } } mitk::FiberBundleX::Pointer fib = dynamic_cast((*it)->GetData()); mitk::FiberBundleX::Pointer newBundle = fib->GetDeepCopy(); QString name((*it)->GetName().c_str()); name += "_copy"; mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); if (parentNode.IsNotNull()) GetDataStorage()->Add(fbNode, parentNode); else GetDataStorage()->Add(fbNode); // copy child fiducials if (m_Controls->m_IncludeFiducials->isChecked()) { mitk::DataStorage::SetOfObjects::ConstPointer derivations = GetDataStorage()->GetDerivations(*it); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = derivations->begin(); it2 != derivations->end(); ++it2 ) { mitk::DataNode::Pointer fiducialNode = *it2; if ( fiducialNode.IsNotNull() && dynamic_cast(fiducialNode->GetData()) ) { mitk::PlanarEllipse::Pointer pe = mitk::PlanarEllipse::New(); pe->DeepCopy(dynamic_cast(fiducialNode->GetData())); mitk::DataNode::Pointer newNode = mitk::DataNode::New(); newNode->SetData(pe); newNode->SetName(fiducialNode->GetName()); newNode->SetBoolProperty("planarfigure.3drendering", true); GetDataStorage()->Add(newNode, fbNode); } } } } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberfoxView::JoinBundles() { if ( m_SelectedBundles.size()<2 ){ QMessageBox::information( NULL, "Warning", "Select at least two fiber bundles!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!"; return; } std::vector::const_iterator it = m_SelectedBundles.begin(); mitk::FiberBundleX::Pointer newBundle = dynamic_cast((*it)->GetData()); QString name(""); name += QString((*it)->GetName().c_str()); ++it; for (; it!=m_SelectedBundles.end(); ++it) { newBundle = newBundle->AddBundle(dynamic_cast((*it)->GetData())); name += "+"+QString((*it)->GetName().c_str()); } mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberfoxView::UpdateGui() { m_Controls->m_FiberBundleLabel->setText("mandatory"); m_Controls->m_GeometryFrame->setEnabled(true); m_Controls->m_GeometryMessage->setVisible(false); m_Controls->m_DiffusionPropsMessage->setVisible(false); m_Controls->m_FiberGenMessage->setVisible(true); m_Controls->m_TransformBundlesButton->setEnabled(false); m_Controls->m_CopyBundlesButton->setEnabled(false); m_Controls->m_GenerateFibersButton->setEnabled(false); m_Controls->m_FlipButton->setEnabled(false); m_Controls->m_CircleButton->setEnabled(false); m_Controls->m_BvalueBox->setEnabled(true); m_Controls->m_NumGradientsBox->setEnabled(true); m_Controls->m_JoinBundlesButton->setEnabled(false); m_Controls->m_AlignOnGrid->setEnabled(false); if (m_SelectedFiducial.IsNotNull()) { m_Controls->m_TransformBundlesButton->setEnabled(true); m_Controls->m_FlipButton->setEnabled(true); m_Controls->m_AlignOnGrid->setEnabled(true); } if (m_SelectedImage.IsNotNull() || !m_SelectedBundles.empty()) { m_Controls->m_TransformBundlesButton->setEnabled(true); m_Controls->m_CircleButton->setEnabled(true); m_Controls->m_FiberGenMessage->setVisible(false); m_Controls->m_AlignOnGrid->setEnabled(true); } if (m_MaskImageNode.IsNotNull() || m_SelectedImage.IsNotNull()) { m_Controls->m_GeometryMessage->setVisible(true); m_Controls->m_GeometryFrame->setEnabled(false); } if (m_SelectedDWI.IsNotNull()) { m_Controls->m_DiffusionPropsMessage->setVisible(true); m_Controls->m_BvalueBox->setEnabled(false); m_Controls->m_NumGradientsBox->setEnabled(false); m_Controls->m_GeometryMessage->setVisible(true); m_Controls->m_GeometryFrame->setEnabled(false); } if (!m_SelectedBundles.empty()) { m_Controls->m_CopyBundlesButton->setEnabled(true); m_Controls->m_GenerateFibersButton->setEnabled(true); m_Controls->m_FiberBundleLabel->setText(m_SelectedBundles.at(0)->GetName().c_str()); if (m_SelectedBundles.size()>1) m_Controls->m_JoinBundlesButton->setEnabled(true); } } void QmitkFiberfoxView::OnSelectionChanged( berry::IWorkbenchPart::Pointer, const QList& nodes ) { m_SelectedBundles2.clear(); m_SelectedImages.clear(); m_SelectedFiducials.clear(); m_SelectedFiducial = NULL; m_SelectedBundles.clear(); m_SelectedImage = NULL; m_SelectedDWI = NULL; m_MaskImageNode = NULL; m_Controls->m_TissueMaskLabel->setText("optional"); // iterate all selected objects, adjust warning visibility for( int i=0; i*>(node->GetData()) ) { m_SelectedDWI = node; m_SelectedImage = node; m_SelectedImages.push_back(node); } else if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { m_SelectedImages.push_back(node); m_SelectedImage = node; bool isbinary = false; node->GetPropertyValue("binary", isbinary); if (isbinary) { m_MaskImageNode = node; m_Controls->m_TissueMaskLabel->setText(m_MaskImageNode->GetName().c_str()); } } else if ( node.IsNotNull() && dynamic_cast(node->GetData()) ) { m_SelectedBundles2.push_back(node); if (m_Controls->m_RealTimeFibers->isChecked()) { m_SelectedBundles.push_back(node); mitk::FiberBundleX::Pointer newFib = dynamic_cast(node->GetData()); if (newFib->GetNumFibers()!=m_Controls->m_FiberDensityBox->value()) GenerateFibers(); } else m_SelectedBundles.push_back(node); } else if ( node.IsNotNull() && dynamic_cast(node->GetData()) ) { m_SelectedFiducials.push_back(node); m_SelectedFiducial = node; m_SelectedBundles.clear(); mitk::DataStorage::SetOfObjects::ConstPointer parents = GetDataStorage()->GetSources(node); for( mitk::DataStorage::SetOfObjects::const_iterator it = parents->begin(); it != parents->end(); ++it ) { mitk::DataNode::Pointer pNode = *it; if ( pNode.IsNotNull() && dynamic_cast(pNode->GetData()) ) m_SelectedBundles.push_back(pNode); } } } UpdateGui(); } void QmitkFiberfoxView::EnableCrosshairNavigation() { MITK_DEBUG << "EnableCrosshairNavigation"; // enable the crosshair navigation if (mitk::ILinkedRenderWindowPart* linkedRenderWindow = dynamic_cast(this->GetRenderWindowPart())) { MITK_DEBUG << "enabling linked navigation"; linkedRenderWindow->EnableLinkedNavigation(true); // linkedRenderWindow->EnableSlicingPlanes(true); } if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::DisableCrosshairNavigation() { MITK_DEBUG << "DisableCrosshairNavigation"; // disable the crosshair navigation during the drawing if (mitk::ILinkedRenderWindowPart* linkedRenderWindow = dynamic_cast(this->GetRenderWindowPart())) { MITK_DEBUG << "disabling linked navigation"; linkedRenderWindow->EnableLinkedNavigation(false); // linkedRenderWindow->EnableSlicingPlanes(false); } } void QmitkFiberfoxView::NodeRemoved(const mitk::DataNode* node) { mitk::DataNode* nonConstNode = const_cast(node); std::map::iterator it = m_DataNodeToPlanarFigureData.find(nonConstNode); if (dynamic_cast(node->GetData())) { m_SelectedBundles.clear(); m_SelectedBundles2.clear(); } else if (dynamic_cast(node->GetData())) m_SelectedImages.clear(); if( it != m_DataNodeToPlanarFigureData.end() ) { QmitkPlanarFigureData& data = it->second; // remove observers data.m_Figure->RemoveObserver( data.m_EndPlacementObserverTag ); data.m_Figure->RemoveObserver( data.m_SelectObserverTag ); data.m_Figure->RemoveObserver( data.m_StartInteractionObserverTag ); data.m_Figure->RemoveObserver( data.m_EndInteractionObserverTag ); m_DataNodeToPlanarFigureData.erase( it ); } } void QmitkFiberfoxView::NodeAdded( const mitk::DataNode* node ) { // add observer for selection in renderwindow mitk::PlanarFigure* figure = dynamic_cast(node->GetData()); bool isPositionMarker (false); node->GetBoolProperty("isContourMarker", isPositionMarker); if( figure && !isPositionMarker ) { MITK_DEBUG << "figure added. will add interactor if needed."; mitk::PlanarFigureInteractor::Pointer figureInteractor = dynamic_cast(node->GetDataInteractor().GetPointer()); mitk::DataNode* nonConstNode = const_cast( node ); if(figureInteractor.IsNull()) { figureInteractor = mitk::PlanarFigureInteractor::New(); us::Module* planarFigureModule = us::ModuleRegistry::GetModule( "MitkPlanarFigure" ); figureInteractor->LoadStateMachine("PlanarFigureInteraction.xml", planarFigureModule ); figureInteractor->SetEventConfig( "PlanarFigureConfig.xml", planarFigureModule ); figureInteractor->SetDataNode( nonConstNode ); } MITK_DEBUG << "will now add observers for planarfigure"; QmitkPlanarFigureData data; data.m_Figure = figure; // // add observer for event when figure has been placed typedef itk::SimpleMemberCommand< QmitkFiberfoxView > SimpleCommandType; // SimpleCommandType::Pointer initializationCommand = SimpleCommandType::New(); // initializationCommand->SetCallbackFunction( this, &QmitkFiberfoxView::PlanarFigureInitialized ); // data.m_EndPlacementObserverTag = figure->AddObserver( mitk::EndPlacementPlanarFigureEvent(), initializationCommand ); // add observer for event when figure is picked (selected) typedef itk::MemberCommand< QmitkFiberfoxView > MemberCommandType; MemberCommandType::Pointer selectCommand = MemberCommandType::New(); selectCommand->SetCallbackFunction( this, &QmitkFiberfoxView::PlanarFigureSelected ); data.m_SelectObserverTag = figure->AddObserver( mitk::SelectPlanarFigureEvent(), selectCommand ); // add observer for event when interaction with figure starts SimpleCommandType::Pointer startInteractionCommand = SimpleCommandType::New(); startInteractionCommand->SetCallbackFunction( this, &QmitkFiberfoxView::DisableCrosshairNavigation); data.m_StartInteractionObserverTag = figure->AddObserver( mitk::StartInteractionPlanarFigureEvent(), startInteractionCommand ); // add observer for event when interaction with figure starts SimpleCommandType::Pointer endInteractionCommand = SimpleCommandType::New(); endInteractionCommand->SetCallbackFunction( this, &QmitkFiberfoxView::EnableCrosshairNavigation); data.m_EndInteractionObserverTag = figure->AddObserver( mitk::EndInteractionPlanarFigureEvent(), endInteractionCommand ); m_DataNodeToPlanarFigureData[nonConstNode] = data; } } void QmitkFiberfoxView::PlanarFigureSelected( itk::Object* object, const itk::EventObject& ) { mitk::TNodePredicateDataType::Pointer isPf = mitk::TNodePredicateDataType::New(); mitk::DataStorage::SetOfObjects::ConstPointer allPfs = this->GetDataStorage()->GetSubset( isPf ); for ( mitk::DataStorage::SetOfObjects::const_iterator it = allPfs->begin(); it!=allPfs->end(); ++it) { mitk::DataNode* node = *it; if( node->GetData() == object ) { node->SetSelected(true); m_SelectedFiducial = node; } else node->SetSelected(false); } UpdateGui(); this->RequestRenderWindowUpdate(); } void QmitkFiberfoxView::SetFocus() { m_Controls->m_CircleButton->setFocus(); } void QmitkFiberfoxView::SetOutputPath() { // SELECT FOLDER DIALOG string outputPath = QFileDialog::getExistingDirectory(NULL, "Save images to...", QString(outputPath.c_str())).toStdString(); if (outputPath.empty()) m_Controls->m_SavePathEdit->setText("-"); else { outputPath += "/"; m_Controls->m_SavePathEdit->setText(QString(outputPath.c_str())); } }