diff --git a/Modules/DiffusionImaging/DiffusionCmdApps/Fiberfox/Fiberfox.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Fiberfox/Fiberfox.cpp index 04e5efbc13..63015abb4f 100755 --- a/Modules/DiffusionImaging/DiffusionCmdApps/Fiberfox/Fiberfox.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Fiberfox/Fiberfox.cpp @@ -1,273 +1,273 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include #include #include using namespace mitk; /*! * \brief Command line interface to Fiberfox. * Simulate a diffusion-weighted image from a tractogram using the specified parameter file. */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Fiberfox"); parser.setCategory("Diffusion Simulation Tools"); parser.setContributor("MIC"); parser.setDescription("Command line interface to Fiberfox." " Simulate a diffusion-weighted image from a tractogram using the specified parameter file."); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "o", mitkCommandLineParser::OutputFile, "Output root:", "output root", us::Any(), false); parser.addArgument("", "i", mitkCommandLineParser::String, "Input:", "Input tractogram or diffusion-weighted image.", us::Any(), false); parser.addArgument("parameters", "p", mitkCommandLineParser::InputFile, "Parameter file:", "fiberfox parameter file (.ffp)", us::Any(), false); parser.addArgument("template", "t", mitkCommandLineParser::String, "Template image:", "Use parameters of the template diffusion-weighted image.", us::Any()); parser.addArgument("verbose", "v", mitkCommandLineParser::Bool, "Output additional images:", "output volume fraction images etc.", us::Any()); parser.addArgument("dont_apply_direction_matrix", "", mitkCommandLineParser::Bool, "Don't apply direction matrix:", "Don't rotate gradients by image direction matrix.", us::Any()); parser.addArgument("fix_seed", "", mitkCommandLineParser::Bool, "Use fix random seed:", "Always use same sequence of random numbers.", us::Any()); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) { return EXIT_FAILURE; } std::string outName = us::any_cast(parsedArgs["o"]); std::string paramName = us::any_cast(parsedArgs["parameters"]); std::string input=""; if (parsedArgs.count("i")) input = us::any_cast(parsedArgs["i"]); bool fix_seed = false; if (parsedArgs.count("fix_seed")) fix_seed = us::any_cast(parsedArgs["fix_seed"]); bool verbose = false; if (parsedArgs.count("verbose")) verbose = us::any_cast(parsedArgs["verbose"]); bool apply_direction_matrix = true; if (parsedArgs.count("dont_apply_direction_matrix")) apply_direction_matrix = false; FiberfoxParameters parameters; - parameters.LoadParameters(paramName); + parameters.LoadParameters(paramName, fix_seed); // Test if /path/dir is an existing directory: std::string file_extension = ""; if( itksys::SystemTools::FileIsDirectory( outName ) ) { while( *(--(outName.cend())) == '/') { outName.pop_back(); } outName = outName + '/'; parameters.m_Misc.m_OutputPath = outName; outName = outName + parameters.m_Misc.m_OutputPrefix; // using default m_OutputPrefix as initialized. } else { // outName is NOT an existing directory, so we need to remove all trailing slashes: while( *(--(outName.cend())) == '/') { outName.pop_back(); } // now split up the given outName into directory and (prefix of) filename: if( ! itksys::SystemTools::GetFilenamePath( outName ).empty() && itksys::SystemTools::FileIsDirectory(itksys::SystemTools::GetFilenamePath( outName ) ) ) { parameters.m_Misc.m_OutputPath = itksys::SystemTools::GetFilenamePath( outName ) + '/'; } else { parameters.m_Misc.m_OutputPath = itksys::SystemTools::GetCurrentWorkingDirectory() + '/'; } file_extension = itksys::SystemTools::GetFilenameExtension(outName); if( ! itksys::SystemTools::GetFilenameWithoutExtension( outName ).empty() ) { parameters.m_Misc.m_OutputPrefix = itksys::SystemTools::GetFilenameWithoutExtension( outName ); } else { parameters.m_Misc.m_OutputPrefix = "fiberfox"; } outName = parameters.m_Misc.m_OutputPath + parameters.m_Misc.m_OutputPrefix; } // check if log file already exists and avoid overwriting existing files: std::string NameTest = outName; int c = 0; while( itksys::SystemTools::FileExists( outName + ".log" ) && c <= std::numeric_limits::max() ) { outName = NameTest + "_" + boost::lexical_cast(c); ++c; } mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images", "Fiberbundles"}, {}); mitk::BaseData::Pointer inputData = mitk::IOUtil::Load(input, &functor)[0]; itk::TractsToDWIImageFilter< short >::Pointer tractsToDwiFilter = itk::TractsToDWIImageFilter< short >::New(); if ( dynamic_cast(inputData.GetPointer()) ) // simulate dataset from fibers { tractsToDwiFilter->SetFiberBundle(dynamic_cast(inputData.GetPointer())); if (parsedArgs.count("template")) { MITK_INFO << "Loading template image"; typedef itk::VectorImage< short, 3 > ItkDwiType; typedef itk::Image< short, 3 > ItkImageType; mitk::BaseData::Pointer templateData = mitk::IOUtil::Load(us::any_cast(parsedArgs["template"]), &functor)[0]; mitk::Image::Pointer template_image = dynamic_cast(templateData.GetPointer()); if (mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(template_image)) { ItkDwiType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::GetItkVectorImage(template_image); parameters.m_SignalGen.m_ImageRegion = itkVectorImagePointer->GetLargestPossibleRegion(); parameters.m_SignalGen.m_ImageSpacing = itkVectorImagePointer->GetSpacing(); parameters.m_SignalGen.m_ImageOrigin = itkVectorImagePointer->GetOrigin(); parameters.m_SignalGen.m_ImageDirection = itkVectorImagePointer->GetDirection(); parameters.SetBvalue(mitk::DiffusionPropertyHelper::GetReferenceBValue(template_image)); parameters.SetGradienDirections(mitk::DiffusionPropertyHelper::GetOriginalGradientContainer(template_image)); } else { ItkImageType::Pointer itkImagePointer = ItkImageType::New(); mitk::CastToItkImage(template_image, itkImagePointer); parameters.m_SignalGen.m_ImageRegion = itkImagePointer->GetLargestPossibleRegion(); parameters.m_SignalGen.m_ImageSpacing = itkImagePointer->GetSpacing(); parameters.m_SignalGen.m_ImageOrigin = itkImagePointer->GetOrigin(); parameters.m_SignalGen.m_ImageDirection = itkImagePointer->GetDirection(); } } } else if ( dynamic_cast(inputData.GetPointer()) ) // add artifacts to existing image { typedef itk::VectorImage< short, 3 > ItkDwiType; mitk::Image::Pointer diffImg = dynamic_cast(inputData.GetPointer()); ItkDwiType::Pointer itkVectorImagePointer = ItkDwiType::New(); mitk::CastToItkImage(diffImg, itkVectorImagePointer); parameters.m_SignalGen.m_SignalScale = 1; parameters.m_SignalGen.m_ImageRegion = itkVectorImagePointer->GetLargestPossibleRegion(); parameters.m_SignalGen.m_ImageSpacing = itkVectorImagePointer->GetSpacing(); parameters.m_SignalGen.m_ImageOrigin = itkVectorImagePointer->GetOrigin(); parameters.m_SignalGen.m_ImageDirection = itkVectorImagePointer->GetDirection(); parameters.SetBvalue(mitk::DiffusionPropertyHelper::GetReferenceBValue(diffImg)); parameters.SetGradienDirections(mitk::DiffusionPropertyHelper::GetOriginalGradientContainer(diffImg)); tractsToDwiFilter->SetInputImage(itkVectorImagePointer); } if (verbose) { MITK_DEBUG << outName << ".ffp"; parameters.SaveParameters(outName+".ffp"); } if (apply_direction_matrix) { MITK_INFO << "Applying direction matrix to gradient directions."; parameters.ApplyDirectionMatrix(); } tractsToDwiFilter->SetParameters(parameters); tractsToDwiFilter->SetUseConstantRandSeed(fix_seed); tractsToDwiFilter->Update(); mitk::Image::Pointer image = mitk::GrabItkImageMemory(tractsToDwiFilter->GetOutput()); if (apply_direction_matrix) mitk::DiffusionPropertyHelper::SetGradientContainer(image, parameters.m_SignalGen.GetItkGradientContainer()); else mitk::DiffusionPropertyHelper::SetOriginalGradientContainer(image, parameters.m_SignalGen.GetItkGradientContainer()); mitk::DiffusionPropertyHelper::SetReferenceBValue(image, parameters.m_SignalGen.GetBvalue()); mitk::DiffusionPropertyHelper::InitializeImage(image); if (file_extension=="") mitk::IOUtil::Save(image, "DWI_NIFTI", outName+".nii.gz"); else if (file_extension==".nii" || file_extension==".nii.gz") mitk::IOUtil::Save(image, "DWI_NIFTI", outName+file_extension); else mitk::IOUtil::Save(image, outName+file_extension); if (verbose) { std::vector< itk::TractsToDWIImageFilter< short >::ItkDoubleImgType::Pointer > volumeFractions = tractsToDwiFilter->GetVolumeFractions(); for (unsigned int k=0; kInitializeByItk(volumeFractions.at(k).GetPointer()); image->SetVolume(volumeFractions.at(k)->GetBufferPointer()); mitk::IOUtil::Save(image, outName+"_Compartment"+boost::lexical_cast(k+1)+".nii.gz"); } if (tractsToDwiFilter->GetPhaseImage().IsNotNull()) { mitk::Image::Pointer image = mitk::Image::New(); itk::TractsToDWIImageFilter< short >::DoubleDwiType::Pointer itkPhase = tractsToDwiFilter->GetPhaseImage(); image = mitk::GrabItkImageMemory( itkPhase.GetPointer() ); mitk::IOUtil::Save(image, outName+"_Phase.nii.gz"); } if (tractsToDwiFilter->GetKspaceImage().IsNotNull()) { mitk::Image::Pointer image = mitk::Image::New(); itk::TractsToDWIImageFilter< short >::DoubleDwiType::Pointer itkImage = tractsToDwiFilter->GetKspaceImage(); image = mitk::GrabItkImageMemory( itkImage.GetPointer() ); mitk::IOUtil::Save(image, outName+"_kSpace.nii.gz"); } int c = 1; std::vector< itk::TractsToDWIImageFilter< short >::DoubleDwiType::Pointer > output_real = tractsToDwiFilter->GetOutputImagesReal(); for (auto real : output_real) { mitk::Image::Pointer image = mitk::Image::New(); image->InitializeByItk(real.GetPointer()); image->SetVolume(real->GetBufferPointer()); mitk::IOUtil::Save(image, outName+"_Coil-"+boost::lexical_cast(c)+"-real.nii.gz"); ++c; } c = 1; std::vector< itk::TractsToDWIImageFilter< short >::DoubleDwiType::Pointer > output_imag = tractsToDwiFilter->GetOutputImagesImag(); for (auto imag : output_imag) { mitk::Image::Pointer image = mitk::Image::New(); image->InitializeByItk(imag.GetPointer()); image->SetVolume(imag->GetBufferPointer()); mitk::IOUtil::Save(image, outName+"_Coil-"+boost::lexical_cast(c)+"-imag.nii.gz"); ++c; } } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp index cebfb6ebf6..24949e8598 100755 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp @@ -1,1770 +1,1770 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "itkTractsToDWIImageFilter.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include namespace itk { template< class PixelType > TractsToDWIImageFilter< PixelType >::TractsToDWIImageFilter() : m_StatusText("") , m_UseConstantRandSeed(false) , m_RandGen(itk::Statistics::MersenneTwisterRandomVariateGenerator::New()) { m_RandGen->SetSeed(); m_DoubleInterpolator = itk::LinearInterpolateImageFunction< ItkDoubleImgType, float >::New(); m_NullDir.Fill(0); } template< class PixelType > TractsToDWIImageFilter< PixelType >::~TractsToDWIImageFilter() { } template< class PixelType > TractsToDWIImageFilter< PixelType >::DoubleDwiType::Pointer TractsToDWIImageFilter< PixelType >:: SimulateKspaceAcquisition( std::vector< DoubleDwiType::Pointer >& compartment_images ) { unsigned int numFiberCompartments = m_Parameters.m_FiberModelList.size(); // create slice object ImageRegion<2> sliceRegion; sliceRegion.SetSize(0, m_WorkingImageRegion.GetSize()[0]); sliceRegion.SetSize(1, m_WorkingImageRegion.GetSize()[1]); Vector< double, 2 > sliceSpacing; sliceSpacing[0] = m_WorkingSpacing[0]; sliceSpacing[1] = m_WorkingSpacing[1]; DoubleDwiType::PixelType nullPix; nullPix.SetSize(compartment_images.at(0)->GetVectorLength()); nullPix.Fill(0.0); auto magnitudeDwiImage = DoubleDwiType::New(); magnitudeDwiImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); magnitudeDwiImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin ); magnitudeDwiImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); magnitudeDwiImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); magnitudeDwiImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); magnitudeDwiImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); magnitudeDwiImage->SetVectorLength( compartment_images.at(0)->GetVectorLength() ); magnitudeDwiImage->Allocate(); magnitudeDwiImage->FillBuffer(nullPix); m_PhaseImage = DoubleDwiType::New(); m_PhaseImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); m_PhaseImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin ); m_PhaseImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); m_PhaseImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_PhaseImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_PhaseImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_PhaseImage->SetVectorLength( compartment_images.at(0)->GetVectorLength() ); m_PhaseImage->Allocate(); m_PhaseImage->FillBuffer(nullPix); m_KspaceImage = DoubleDwiType::New(); m_KspaceImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); m_KspaceImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin ); m_KspaceImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); m_KspaceImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_KspaceImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_KspaceImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_KspaceImage->SetVectorLength( m_Parameters.m_SignalGen.m_NumberOfCoils ); m_KspaceImage->Allocate(); m_KspaceImage->FillBuffer(nullPix); std::list< unsigned int > spikeVolume; if (m_Parameters.m_Misc.m_DoAddSpikes) for (unsigned int i=0; iGetIntegerVariate()%(compartment_images.at(0)->GetVectorLength())); // calculate coil positions double a = m_Parameters.m_SignalGen.m_ImageRegion.GetSize(0)*m_Parameters.m_SignalGen.m_ImageSpacing[0]; double b = m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1)*m_Parameters.m_SignalGen.m_ImageSpacing[1]; double c = m_Parameters.m_SignalGen.m_ImageRegion.GetSize(2)*m_Parameters.m_SignalGen.m_ImageSpacing[2]; double diagonal = sqrt(a*a+b*b)/1000; // image diagonal in m m_CoilPointset = mitk::PointSet::New(); std::vector< itk::Vector > coilPositions; itk::Vector pos; pos.Fill(0.0); pos[1] = -diagonal/2; itk::Vector center; center[0] = a/2-m_Parameters.m_SignalGen.m_ImageSpacing[0]/2; center[1] = b/2-m_Parameters.m_SignalGen.m_ImageSpacing[2]/2; center[2] = c/2-m_Parameters.m_SignalGen.m_ImageSpacing[1]/2; for (int c=0; cInsertPoint(c, pos*1000 + m_Parameters.m_SignalGen.m_ImageOrigin.GetVectorFromOrigin() + center ); double rz = 360.0/m_Parameters.m_SignalGen.m_NumberOfCoils * itk::Math::pi/180; vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); rotZ[0][0] = cos(rz); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(rz); rotZ[1][0] = -rotZ[0][1]; pos.SetVnlVector(rotZ*pos.GetVnlVector()); } auto max_threads = omp_get_max_threads(); int out_threads = Math::ceil(std::sqrt(max_threads)); int in_threads = Math::floor(std::sqrt(max_threads)); - PrintToLog("Parallel volumes: " + boost::lexical_cast(out_threads), false, true, false); - PrintToLog("Threads per slice: " + boost::lexical_cast(in_threads), false, true, false); + PrintToLog("Parallel volumes: " + boost::lexical_cast(out_threads), false, true, true); + PrintToLog("Threads per slice: " + boost::lexical_cast(in_threads), false, true, true); PrintToLog("0% 10 20 30 40 50 60 70 80 90 100%", false, true, false); PrintToLog("|----|----|----|----|----|----|----|----|----|----|\n*", false, false, false); unsigned long lastTick = 0; boost::progress_display disp(compartment_images.at(0)->GetVectorLength()*compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2)); #pragma omp parallel for num_threads(out_threads) for (int g=0; g(compartment_images.at(0)->GetVectorLength()); g++) { if (this->GetAbortGenerateData()) continue; std::list< unsigned int > spikeSlice; #pragma omp critical { for (auto sv : spikeVolume) if (sv == static_cast(g)) spikeSlice.push_back(m_RandGen->GetIntegerVariate()%compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2)); spikeVolume.remove(g); } for (unsigned int z=0; zGetLargestPossibleRegion().GetSize(2); z++) { std::vector< Float2DImageType::Pointer > compartment_slices; std::vector< float > t2Vector; std::vector< float > t1Vector; for (unsigned int i=0; i* signalModel; if (iSetLargestPossibleRegion( sliceRegion ); slice->SetBufferedRegion( sliceRegion ); slice->SetRequestedRegion( sliceRegion ); slice->SetSpacing(sliceSpacing); slice->Allocate(); slice->FillBuffer(0.0); // extract slice from channel g for (unsigned int y=0; yGetLargestPossibleRegion().GetSize(1); y++) for (unsigned int x=0; xGetLargestPossibleRegion().GetSize(0); x++) { Float2DImageType::IndexType index2D; index2D[0]=x; index2D[1]=y; DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z; slice->SetPixel(index2D, compartment_images.at(i)->GetPixel(index3D)[g]); } compartment_slices.push_back(slice); t2Vector.push_back(signalModel->GetT2()); t1Vector.push_back(signalModel->GetT1()); } int numSpikes = 0; for (auto ss : spikeSlice) if (ss == z) ++numSpikes; spikeSlice.remove(z); int spikeCoil = m_RandGen->GetIntegerVariate()%m_Parameters.m_SignalGen.m_NumberOfCoils; if (this->GetAbortGenerateData()) continue; for (int c=0; c::New(); idft->SetCompartmentImages(compartment_slices); idft->SetT2(t2Vector); idft->SetT1(t1Vector); idft->SetUseConstantRandSeed(m_UseConstantRandSeed); idft->SetParameters(&m_Parameters); idft->SetZ((float)z-(float)( compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2) -compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2)%2 ) / 2.0); idft->SetZidx(z); idft->SetCoilPosition(coilPositions.at(c)); idft->SetFiberBundle(m_FiberBundle); idft->SetTranslation(m_Translations.at(g)); idft->SetRotationMatrix(m_RotationsInv.at(g)); idft->SetDiffusionGradientDirection(m_Parameters.m_SignalGen.GetGradientDirection(g)); if (c==spikeCoil) idft->SetSpikesPerSlice(numSpikes); idft->SetNumberOfThreads(in_threads); idft->Update(); #pragma omp critical if (c==spikeCoil && numSpikes>0) { m_SpikeLog += "Volume " + boost::lexical_cast(g) + " Coil " + boost::lexical_cast(c) + "\n"; m_SpikeLog += idft->GetSpikeLog(); } Complex2DImageType::Pointer fSlice; fSlice = idft->GetOutput(); // fourier transform slice Complex2DImageType::Pointer newSlice; auto dft = itk::DftImageFilter< Float2DImageType::PixelType >::New(); dft->SetInput(fSlice); dft->SetParameters(m_Parameters); dft->SetNumberOfThreads(in_threads); dft->Update(); newSlice = dft->GetOutput(); // put slice back into channel g for (unsigned int y=0; yGetLargestPossibleRegion().GetSize(1); y++) for (unsigned int x=0; xGetLargestPossibleRegion().GetSize(0); x++) { DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z; Complex2DImageType::IndexType index2D; index2D[0]=x; index2D[1]=y; Complex2DImageType::PixelType cPix = newSlice->GetPixel(index2D); double magn = sqrt(cPix.real()*cPix.real()+cPix.imag()*cPix.imag()); double phase = 0; if (cPix.real()!=0) phase = atan( cPix.imag()/cPix.real() ); DoubleDwiType::PixelType real_pix = m_OutputImagesReal.at(c)->GetPixel(index3D); real_pix[g] = cPix.real(); m_OutputImagesReal.at(c)->SetPixel(index3D, real_pix); DoubleDwiType::PixelType imag_pix = m_OutputImagesImag.at(c)->GetPixel(index3D); imag_pix[g] = cPix.imag(); m_OutputImagesImag.at(c)->SetPixel(index3D, imag_pix); DoubleDwiType::PixelType dwiPix = magnitudeDwiImage->GetPixel(index3D); DoubleDwiType::PixelType phasePix = m_PhaseImage->GetPixel(index3D); if (m_Parameters.m_SignalGen.m_NumberOfCoils>1) { dwiPix[g] += magn*magn; phasePix[g] += phase*phase; } else { dwiPix[g] = magn; phasePix[g] = phase; } //#pragma omp critical { magnitudeDwiImage->SetPixel(index3D, dwiPix); m_PhaseImage->SetPixel(index3D, phasePix); // k-space image if (g==0) { DoubleDwiType::PixelType kspacePix = m_KspaceImage->GetPixel(index3D); kspacePix[c] = idft->GetKSpaceImage()->GetPixel(index2D); m_KspaceImage->SetPixel(index3D, kspacePix); } } } } if (m_Parameters.m_SignalGen.m_NumberOfCoils>1) { for (int y=0; y(magnitudeDwiImage->GetLargestPossibleRegion().GetSize(1)); y++) for (int x=0; x(magnitudeDwiImage->GetLargestPossibleRegion().GetSize(0)); x++) { DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z; DoubleDwiType::PixelType magPix = magnitudeDwiImage->GetPixel(index3D); magPix[g] = sqrt(magPix[g]/m_Parameters.m_SignalGen.m_NumberOfCoils); DoubleDwiType::PixelType phasePix = m_PhaseImage->GetPixel(index3D); phasePix[g] = sqrt(phasePix[g]/m_Parameters.m_SignalGen.m_NumberOfCoils); //#pragma omp critical { magnitudeDwiImage->SetPixel(index3D, magPix); m_PhaseImage->SetPixel(index3D, phasePix); } } } ++disp; unsigned long newTick = 50*disp.count()/disp.expected_count(); for (unsigned long tick = 0; tick<(newTick-lastTick); tick++) PrintToLog("*", false, false, false); lastTick = newTick; } } PrintToLog("\n", false); return magnitudeDwiImage; } template< class PixelType > TractsToDWIImageFilter< PixelType >::ItkDoubleImgType::Pointer TractsToDWIImageFilter< PixelType >:: NormalizeInsideMask(ItkDoubleImgType::Pointer image) { double max = itk::NumericTraits< double >::min(); double min = itk::NumericTraits< double >::max(); itk::ImageRegionIterator< ItkDoubleImgType > it(image, image->GetLargestPossibleRegion()); while(!it.IsAtEnd()) { if (m_Parameters.m_SignalGen.m_MaskImage.IsNotNull() && m_Parameters.m_SignalGen.m_MaskImage->GetPixel(it.GetIndex())<=0) { it.Set(0.0); ++it; continue; } if (it.Get()>max) max = it.Get(); if (it.Get()::New(); scaler->SetInput(image); scaler->SetShift(-min); scaler->SetScale(1.0/(max-min)); scaler->Update(); return scaler->GetOutput(); } template< class PixelType > void TractsToDWIImageFilter< PixelType >::CheckVolumeFractionImages() { m_UseRelativeNonFiberVolumeFractions = false; // check for fiber volume fraction maps unsigned int fibVolImages = 0; for (std::size_t i=0; iGetVolumeFractionImage().IsNotNull()) { PrintToLog("Using volume fraction map for fiber compartment " + boost::lexical_cast(i+1), false); fibVolImages++; } } // check for non-fiber volume fraction maps unsigned int nonfibVolImages = 0; for (std::size_t i=0; iGetVolumeFractionImage().IsNotNull()) { PrintToLog("Using volume fraction map for non-fiber compartment " + boost::lexical_cast(i+1), false); nonfibVolImages++; } } // not all fiber compartments are using volume fraction maps // --> non-fiber volume fractions are assumed to be relative to the // non-fiber volume and not absolute voxel-volume fractions. // this means if two non-fiber compartments are used but only one of them // has an associated volume fraction map, the repesctive other volume fraction map // can be determined as inverse (1-val) of the present volume fraction map- if ( fibVolImages::New(); inverter->SetMaximum(1.0); if ( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage().IsNull() && m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage().IsNotNull() ) { // m_Parameters.m_NonFiberModelList[1]->SetVolumeFractionImage( // NormalizeInsideMask( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage() ) ); inverter->SetInput( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage() ); inverter->Update(); m_Parameters.m_NonFiberModelList[0]->SetVolumeFractionImage(inverter->GetOutput()); } else if ( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage().IsNull() && m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage().IsNotNull() ) { // m_Parameters.m_NonFiberModelList[0]->SetVolumeFractionImage( // NormalizeInsideMask( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage() ) ); inverter->SetInput( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage() ); inverter->Update(); m_Parameters.m_NonFiberModelList[1]->SetVolumeFractionImage(inverter->GetOutput()); } else { itkExceptionMacro("Something went wrong in automatically calculating the missing non-fiber volume fraction image!" " Did you use two non fiber compartments but only one volume fraction image?" " Then it should work and this error is really strange."); } m_UseRelativeNonFiberVolumeFractions = true; nonfibVolImages++; } // Up to two fiber compartments are allowed without volume fraction maps since the volume fractions can then be determined automatically if (m_Parameters.m_FiberModelList.size()>2 && fibVolImages!=m_Parameters.m_FiberModelList.size()) itkExceptionMacro("More than two fiber compartment selected but no corresponding volume fraction maps set!"); // One non-fiber compartment is allowed without volume fraction map since the volume fraction can then be determined automatically if (m_Parameters.m_NonFiberModelList.size()>1 && nonfibVolImages!=m_Parameters.m_NonFiberModelList.size()) itkExceptionMacro("More than one non-fiber compartment selected but no volume fraction maps set!"); if (fibVolImages0) { PrintToLog("Not all fiber compartments are using an associated volume fraction image.\n" "Assuming non-fiber volume fraction images to contain values relative to the" " remaining non-fiber volume, not absolute values.", false); m_UseRelativeNonFiberVolumeFractions = true; // mitk::LocaleSwitch localeSwitch("C"); // itk::ImageFileWriter::Pointer wr = itk::ImageFileWriter::New(); // wr->SetInput(m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage()); // wr->SetFileName("/local/volumefraction.nrrd"); // wr->Update(); } // initialize the images that store the output volume fraction of each compartment m_VolumeFractions.clear(); for (std::size_t i=0; iSetSpacing( m_WorkingSpacing ); doubleImg->SetOrigin( m_WorkingOrigin ); doubleImg->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); doubleImg->SetLargestPossibleRegion( m_WorkingImageRegion ); doubleImg->SetBufferedRegion( m_WorkingImageRegion ); doubleImg->SetRequestedRegion( m_WorkingImageRegion ); doubleImg->Allocate(); doubleImg->FillBuffer(0); m_VolumeFractions.push_back(doubleImg); } } template< class PixelType > void TractsToDWIImageFilter< PixelType >::InitializeData() { m_Rotations.clear(); m_Translations.clear(); m_MotionLog = ""; m_SpikeLog = ""; // initialize output dwi image m_Parameters.m_SignalGen.m_CroppedRegion = m_Parameters.m_SignalGen.m_ImageRegion; if (m_Parameters.m_Misc.m_DoAddAliasing) m_Parameters.m_SignalGen.m_CroppedRegion.SetSize( 1, m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(1) *m_Parameters.m_SignalGen.m_CroppingFactor); itk::Point shiftedOrigin = m_Parameters.m_SignalGen.m_ImageOrigin; shiftedOrigin[1] += (m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1) -m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(1))*m_Parameters.m_SignalGen.m_ImageSpacing[1]/2; m_OutputImage = OutputImageType::New(); m_OutputImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); m_OutputImage->SetOrigin( shiftedOrigin ); m_OutputImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); m_OutputImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_OutputImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_OutputImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_OutputImage->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() ); m_OutputImage->Allocate(); typename OutputImageType::PixelType temp; temp.SetSize(m_Parameters.m_SignalGen.GetNumVolumes()); temp.Fill(0.0); m_OutputImage->FillBuffer(temp); PrintToLog("Output image spacing: [" + boost::lexical_cast(m_Parameters.m_SignalGen.m_ImageSpacing[0]) + "," + boost::lexical_cast(m_Parameters.m_SignalGen.m_ImageSpacing[1]) + "," + boost::lexical_cast(m_Parameters.m_SignalGen.m_ImageSpacing[2]) + "]", false); PrintToLog("Output image size: [" + boost::lexical_cast(m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(0)) + "," + boost::lexical_cast(m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(1)) + "," + boost::lexical_cast(m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(2)) + "]", false); // images containing real and imaginary part of the dMRI signal for each coil m_OutputImagesReal.clear(); m_OutputImagesImag.clear(); for (int i=0; iSetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); outputImageReal->SetOrigin( shiftedOrigin ); outputImageReal->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); outputImageReal->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); outputImageReal->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); outputImageReal->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); outputImageReal->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() ); outputImageReal->Allocate(); outputImageReal->FillBuffer(temp); m_OutputImagesReal.push_back(outputImageReal); typename DoubleDwiType::Pointer outputImageImag = DoubleDwiType::New(); outputImageImag->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); outputImageImag->SetOrigin( shiftedOrigin ); outputImageImag->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); outputImageImag->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); outputImageImag->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); outputImageImag->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); outputImageImag->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() ); outputImageImag->Allocate(); outputImageImag->FillBuffer(temp); m_OutputImagesImag.push_back(outputImageImag); } // Apply in-plane upsampling for Gibbs ringing artifact double upsampling = 1; if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging) upsampling = 2; m_WorkingSpacing = m_Parameters.m_SignalGen.m_ImageSpacing; m_WorkingSpacing[0] /= upsampling; m_WorkingSpacing[1] /= upsampling; m_WorkingImageRegion = m_Parameters.m_SignalGen.m_ImageRegion; m_WorkingImageRegion.SetSize(0, m_Parameters.m_SignalGen.m_ImageRegion.GetSize()[0]*upsampling); m_WorkingImageRegion.SetSize(1, m_Parameters.m_SignalGen.m_ImageRegion.GetSize()[1]*upsampling); m_WorkingOrigin = m_Parameters.m_SignalGen.m_ImageOrigin; m_WorkingOrigin[0] -= m_Parameters.m_SignalGen.m_ImageSpacing[0]/2; m_WorkingOrigin[0] += m_WorkingSpacing[0]/2; m_WorkingOrigin[1] -= m_Parameters.m_SignalGen.m_ImageSpacing[1]/2; m_WorkingOrigin[1] += m_WorkingSpacing[1]/2; m_WorkingOrigin[2] -= m_Parameters.m_SignalGen.m_ImageSpacing[2]/2; m_WorkingOrigin[2] += m_WorkingSpacing[2]/2; m_VoxelVolume = m_WorkingSpacing[0]*m_WorkingSpacing[1]*m_WorkingSpacing[2]; PrintToLog("Working image spacing: [" + boost::lexical_cast(m_WorkingSpacing[0]) + "," + boost::lexical_cast(m_WorkingSpacing[1]) + "," + boost::lexical_cast(m_WorkingSpacing[2]) + "]", false); PrintToLog("Working image size: [" + boost::lexical_cast(m_WorkingImageRegion.GetSize(0)) + "," + boost::lexical_cast(m_WorkingImageRegion.GetSize(1)) + "," + boost::lexical_cast(m_WorkingImageRegion.GetSize(2)) + "]", false); // generate double images to store the individual compartment signals m_CompartmentImages.clear(); int numFiberCompartments = m_Parameters.m_FiberModelList.size(); int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size(); for (int i=0; iSetSpacing( m_WorkingSpacing ); doubleDwi->SetOrigin( m_WorkingOrigin ); doubleDwi->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); doubleDwi->SetLargestPossibleRegion( m_WorkingImageRegion ); doubleDwi->SetBufferedRegion( m_WorkingImageRegion ); doubleDwi->SetRequestedRegion( m_WorkingImageRegion ); doubleDwi->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() ); doubleDwi->Allocate(); DoubleDwiType::PixelType pix; pix.SetSize(m_Parameters.m_SignalGen.GetNumVolumes()); pix.Fill(0.0); doubleDwi->FillBuffer(pix); m_CompartmentImages.push_back(doubleDwi); } if (m_FiberBundle.IsNull() && m_InputImage.IsNotNull()) { m_CompartmentImages.clear(); m_Parameters.m_SignalGen.m_DoAddMotion = false; m_Parameters.m_SignalGen.m_DoSimulateRelaxation = false; PrintToLog("Simulating acquisition for input diffusion-weighted image.", false); auto caster = itk::CastImageFilter< OutputImageType, DoubleDwiType >::New(); caster->SetInput(m_InputImage); caster->Update(); if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging) { PrintToLog("Upsampling input diffusion-weighted image for Gibbs ringing simulation.", false); auto resampler = itk::ResampleDwiImageFilter< double >::New(); resampler->SetInput(caster->GetOutput()); itk::Vector< double, 3 > samplingFactor; samplingFactor[0] = upsampling; samplingFactor[1] = upsampling; samplingFactor[2] = 1; resampler->SetSamplingFactor(samplingFactor); resampler->SetInterpolation(itk::ResampleDwiImageFilter< double >::Interpolate_WindowedSinc); resampler->Update(); m_CompartmentImages.push_back(resampler->GetOutput()); } else m_CompartmentImages.push_back(caster->GetOutput()); VectorType translation; translation.Fill(0.0); MatrixType rotation; rotation.SetIdentity(); for (unsigned int g=0; gGetLargestPossibleRegion()!=m_WorkingImageRegion) { PrintToLog("Resampling tissue mask", false); // rescale mask image (otherwise there are problems with the resampling) auto rescaler = itk::RescaleIntensityImageFilter::New(); rescaler->SetInput(0,m_Parameters.m_SignalGen.m_MaskImage); rescaler->SetOutputMaximum(100); rescaler->SetOutputMinimum(0); rescaler->Update(); // resample mask image auto resampler = itk::ResampleImageFilter::New(); resampler->SetInput(rescaler->GetOutput()); resampler->SetSize(m_WorkingImageRegion.GetSize()); resampler->SetOutputSpacing(m_WorkingSpacing); resampler->SetOutputOrigin(m_WorkingOrigin); resampler->SetOutputDirection(m_Parameters.m_SignalGen.m_ImageDirection); resampler->SetOutputStartIndex ( m_WorkingImageRegion.GetIndex() ); auto nn_interpolator = itk::NearestNeighborInterpolateImageFunction::New(); resampler->SetInterpolator(nn_interpolator); resampler->Update(); m_Parameters.m_SignalGen.m_MaskImage = resampler->GetOutput(); } // resample frequency map if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull() && m_Parameters.m_SignalGen.m_FrequencyMap->GetLargestPossibleRegion()!=m_WorkingImageRegion) { PrintToLog("Resampling frequency map", false); auto resampler = itk::ResampleImageFilter::New(); resampler->SetInput(m_Parameters.m_SignalGen.m_FrequencyMap); resampler->SetSize(m_WorkingImageRegion.GetSize()); resampler->SetOutputSpacing(m_WorkingSpacing); resampler->SetOutputOrigin(m_WorkingOrigin); resampler->SetOutputDirection(m_Parameters.m_SignalGen.m_ImageDirection); resampler->SetOutputStartIndex ( m_WorkingImageRegion.GetIndex() ); auto nn_interpolator = itk::NearestNeighborInterpolateImageFunction::New(); resampler->SetInterpolator(nn_interpolator); resampler->Update(); m_Parameters.m_SignalGen.m_FrequencyMap = resampler->GetOutput(); } m_MaskImageSet = true; if (m_Parameters.m_SignalGen.m_MaskImage.IsNull()) { // no input tissue mask is set -> create default PrintToLog("No tissue mask set", false); m_Parameters.m_SignalGen.m_MaskImage = ItkUcharImgType::New(); m_Parameters.m_SignalGen.m_MaskImage->SetSpacing( m_WorkingSpacing ); m_Parameters.m_SignalGen.m_MaskImage->SetOrigin( m_WorkingOrigin ); m_Parameters.m_SignalGen.m_MaskImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); m_Parameters.m_SignalGen.m_MaskImage->SetLargestPossibleRegion( m_WorkingImageRegion ); m_Parameters.m_SignalGen.m_MaskImage->SetBufferedRegion( m_WorkingImageRegion ); m_Parameters.m_SignalGen.m_MaskImage->SetRequestedRegion( m_WorkingImageRegion ); m_Parameters.m_SignalGen.m_MaskImage->Allocate(); m_Parameters.m_SignalGen.m_MaskImage->FillBuffer(100); m_MaskImageSet = false; } else { PrintToLog("Using tissue mask", false); } if (m_Parameters.m_SignalGen.m_DoAddMotion) { if (m_Parameters.m_SignalGen.m_DoRandomizeMotion) { PrintToLog("Random motion artifacts:", false); PrintToLog("Maximum rotation: +/-" + boost::lexical_cast(m_Parameters.m_SignalGen.m_Rotation) + "°", false); PrintToLog("Maximum translation: +/-" + boost::lexical_cast(m_Parameters.m_SignalGen.m_Translation) + "mm", false); } else { PrintToLog("Linear motion artifacts:", false); PrintToLog("Maximum rotation: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_Rotation) + "°", false); PrintToLog("Maximum translation: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_Translation) + "mm", false); } } if ( m_Parameters.m_SignalGen.m_MotionVolumes.empty() ) { // no motion in first volume m_Parameters.m_SignalGen.m_MotionVolumes.push_back(false); // motion in all other volumes while ( m_Parameters.m_SignalGen.m_MotionVolumes.size() < m_Parameters.m_SignalGen.GetNumVolumes() ) { m_Parameters.m_SignalGen.m_MotionVolumes.push_back(true); } } // we need to know for every volume if there is motion. if this information is missing, then set corresponding fal to false while ( m_Parameters.m_SignalGen.m_MotionVolumes.size()::New(); duplicator->SetInputImage(m_Parameters.m_SignalGen.m_MaskImage); duplicator->Update(); m_TransformedMaskImage = duplicator->GetOutput(); // second upsampling needed for motion artifacts ImageRegion<3> upsampledImageRegion = m_WorkingImageRegion; VectorType upsampledSpacing = m_WorkingSpacing; upsampledSpacing[0] /= 4; upsampledSpacing[1] /= 4; upsampledSpacing[2] /= 4; upsampledImageRegion.SetSize(0, m_WorkingImageRegion.GetSize()[0]*4); upsampledImageRegion.SetSize(1, m_WorkingImageRegion.GetSize()[1]*4); upsampledImageRegion.SetSize(2, m_WorkingImageRegion.GetSize()[2]*4); itk::Point upsampledOrigin = m_WorkingOrigin; upsampledOrigin[0] -= m_WorkingSpacing[0]/2; upsampledOrigin[0] += upsampledSpacing[0]/2; upsampledOrigin[1] -= m_WorkingSpacing[1]/2; upsampledOrigin[1] += upsampledSpacing[1]/2; upsampledOrigin[2] -= m_WorkingSpacing[2]/2; upsampledOrigin[2] += upsampledSpacing[2]/2; m_UpsampledMaskImage = ItkUcharImgType::New(); auto upsampler = itk::ResampleImageFilter::New(); upsampler->SetInput(m_Parameters.m_SignalGen.m_MaskImage); upsampler->SetOutputParametersFromImage(m_Parameters.m_SignalGen.m_MaskImage); upsampler->SetSize(upsampledImageRegion.GetSize()); upsampler->SetOutputSpacing(upsampledSpacing); upsampler->SetOutputOrigin(upsampledOrigin); auto nn_interpolator = itk::NearestNeighborInterpolateImageFunction::New(); upsampler->SetInterpolator(nn_interpolator); upsampler->Update(); m_UpsampledMaskImage = upsampler->GetOutput(); } template< class PixelType > void TractsToDWIImageFilter< PixelType >::InitializeFiberData() { m_mmRadius = m_Parameters.m_SignalGen.m_AxonRadius/1000; auto caster = itk::CastImageFilter< itk::Image, itk::Image >::New(); caster->SetInput(m_TransformedMaskImage); caster->Update(); vtkSmartPointer weights = m_FiberBundle->GetFiberWeights(); float mean_weight = 0; for (int i=0; iGetSize(); i++) mean_weight += weights->GetValue(i); mean_weight /= weights->GetSize(); if (mean_weight>0.000001) for (int i=0; iGetSize(); i++) m_FiberBundle->SetFiberWeight(i, weights->GetValue(i)/mean_weight); else PrintToLog("\nWarning: streamlines have VERY low weights. Average weight: " + boost::lexical_cast(mean_weight) + ". Possible source of calculation errors.", false, true, true); auto density_calculator = itk::TractDensityImageFilter< itk::Image >::New(); density_calculator->SetFiberBundle(m_FiberBundle); density_calculator->SetInputImage(caster->GetOutput()); density_calculator->SetBinaryOutput(false); density_calculator->SetUseImageGeometry(true); density_calculator->SetOutputAbsoluteValues(true); density_calculator->Update(); double max_density = density_calculator->GetMaxDensity(); double voxel_volume = m_WorkingSpacing[0]*m_WorkingSpacing[1]*m_WorkingSpacing[2]; if (m_mmRadius>0) { std::stringstream stream; stream << std::fixed << setprecision(2) << itk::Math::pi*m_mmRadius*m_mmRadius*max_density; std::string s = stream.str(); PrintToLog("\nMax. fiber volume: " + s + "mm².", false, true, true); { double full_radius = 1000*std::sqrt(voxel_volume/(max_density*itk::Math::pi)); std::stringstream stream; stream << std::fixed << setprecision(2) << full_radius; std::string s = stream.str(); PrintToLog("\nA full fiber voxel corresponds to a fiber radius of ~" + s + "µm, given the current fiber configuration.", false, true, true); } } else { m_mmRadius = std::sqrt(voxel_volume/(max_density*itk::Math::pi)); std::stringstream stream; stream << std::fixed << setprecision(2) << m_mmRadius*1000; std::string s = stream.str(); PrintToLog("\nSetting fiber radius to " + s + "µm to obtain full voxel.", false, true, true); } // a second fiber bundle is needed to store the transformed version of the m_FiberBundleWorkingCopy m_FiberBundleTransformed = m_FiberBundle->GetDeepCopy(); } template< class PixelType > bool TractsToDWIImageFilter< PixelType >::PrepareLogFile() { if(m_Logfile.is_open()) m_Logfile.close(); std::string filePath; std::string fileName; // Get directory name: if (m_Parameters.m_Misc.m_OutputPath.size() > 0) { filePath = m_Parameters.m_Misc.m_OutputPath; if( *(--(filePath.cend())) != '/') { filePath.push_back('/'); } } else { filePath = mitk::IOUtil::GetTempPath() + '/'; } // check if directory exists, else use /tmp/: if( itksys::SystemTools::FileIsDirectory( filePath ) ) { while( *(--(filePath.cend())) == '/') { filePath.pop_back(); } filePath = filePath + '/'; } else { filePath = mitk::IOUtil::GetTempPath() + '/'; } // Get file name: if( ! m_Parameters.m_Misc.m_ResultNode->GetName().empty() ) { fileName = m_Parameters.m_Misc.m_ResultNode->GetName(); } else { fileName = ""; } if( ! m_Parameters.m_Misc.m_OutputPrefix.empty() ) { fileName = m_Parameters.m_Misc.m_OutputPrefix + fileName; } else { fileName = "fiberfox"; } // check if file already exists and DO NOT overwrite existing files: std::string NameTest = fileName; int c = 0; while( itksys::SystemTools::FileExists( filePath + '/' + fileName + ".log" ) && c <= std::numeric_limits::max() ) { fileName = NameTest + "_" + boost::lexical_cast(c); ++c; } try { m_Logfile.open( ( filePath + '/' + fileName + ".log" ).c_str() ); } catch (const std::ios_base::failure &fail) { MITK_ERROR << "itkTractsToDWIImageFilter.cpp: Exception " << fail.what() << " while trying to open file" << filePath << '/' << fileName << ".log"; return false; } if ( m_Logfile.is_open() ) { PrintToLog( "Logfile: " + filePath + '/' + fileName + ".log", false ); return true; } else { m_StatusText += "Logfile could not be opened!\n"; MITK_ERROR << "itkTractsToDWIImageFilter.cpp: Logfile could not be opened!"; return false; } } template< class PixelType > void TractsToDWIImageFilter< PixelType >::GenerateData() { PrintToLog("\n**********************************************", false); // prepare logfile if ( ! PrepareLogFile() ) { this->SetAbortGenerateData( true ); return; } PrintToLog("Starting Fiberfox dMRI simulation"); m_TimeProbe.Start(); // check input data if (m_FiberBundle.IsNull() && m_InputImage.IsNull()) itkExceptionMacro("Input fiber bundle and input diffusion-weighted image is nullptr!"); if (m_Parameters.m_FiberModelList.empty() && m_InputImage.IsNull()) itkExceptionMacro("No diffusion model for fiber compartments defined and input diffusion-weighted" " image is nullptr! At least one fiber compartment is necessary to simulate diffusion."); if (m_Parameters.m_NonFiberModelList.empty() && m_InputImage.IsNull()) itkExceptionMacro("No diffusion model for non-fiber compartments defined and input diffusion-weighted" " image is nullptr! At least one non-fiber compartment is necessary to simulate diffusion."); if (m_Parameters.m_SignalGen.m_DoDisablePartialVolume) // no partial volume? remove all but first fiber compartment while (m_Parameters.m_FiberModelList.size()>1) m_Parameters.m_FiberModelList.pop_back(); if (!m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition) // No upsampling of input image needed if no k-space simulation is performed m_Parameters.m_SignalGen.m_DoAddGibbsRinging = false; if (m_UseConstantRandSeed) // always generate the same random numbers? m_RandGen->SetSeed(0); else m_RandGen->SetSeed(); InitializeData(); if ( m_FiberBundle.IsNotNull() ) // if no fiber bundle is found, we directly proceed to the k-space acquisition simulation { CheckVolumeFractionImages(); InitializeFiberData(); int numFiberCompartments = m_Parameters.m_FiberModelList.size(); int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size(); double maxVolume = 0; unsigned long lastTick = 0; int signalModelSeed = m_RandGen->GetIntegerVariate(); PrintToLog("\n", false, false); PrintToLog("Generating " + boost::lexical_cast(numFiberCompartments+numNonFiberCompartments) + "-compartment diffusion-weighted signal."); std::vector< int > bVals = m_Parameters.m_SignalGen.GetBvalues(); PrintToLog("b-values: ", false, false, true); for (auto v : bVals) PrintToLog(boost::lexical_cast(v) + " ", false, false, true); PrintToLog("\nVolumes: " + boost::lexical_cast(m_Parameters.m_SignalGen.GetNumVolumes()), false, true, true); PrintToLog("\n", false, false, true); PrintToLog("\n", false, false, true); unsigned int image_size_x = m_WorkingImageRegion.GetSize(0); unsigned int region_size_y = m_WorkingImageRegion.GetSize(1); unsigned int num_gradients = m_Parameters.m_SignalGen.GetNumVolumes(); int numFibers = m_FiberBundle->GetNumFibers(); boost::progress_display disp(numFibers*num_gradients); if (m_FiberBundle->GetMeanFiberLength()<5.0) omp_set_num_threads(2); PrintToLog("0% 10 20 30 40 50 60 70 80 90 100%", false, true, false); PrintToLog("|----|----|----|----|----|----|----|----|----|----|\n*", false, false, false); for (unsigned int g=0; gSetSeed(signalModelSeed); for (std::size_t i=0; iSetSeed(signalModelSeed); // storing voxel-wise intra-axonal volume in mm³ auto intraAxonalVolumeImage = ItkDoubleImgType::New(); intraAxonalVolumeImage->SetSpacing( m_WorkingSpacing ); intraAxonalVolumeImage->SetOrigin( m_WorkingOrigin ); intraAxonalVolumeImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); intraAxonalVolumeImage->SetLargestPossibleRegion( m_WorkingImageRegion ); intraAxonalVolumeImage->SetBufferedRegion( m_WorkingImageRegion ); intraAxonalVolumeImage->SetRequestedRegion( m_WorkingImageRegion ); intraAxonalVolumeImage->Allocate(); intraAxonalVolumeImage->FillBuffer(0); maxVolume = 0; double* intraAxBuffer = intraAxonalVolumeImage->GetBufferPointer(); if (this->GetAbortGenerateData()) continue; vtkPolyData* fiberPolyData = m_FiberBundleTransformed->GetFiberPolyData(); // generate fiber signal (if there are any fiber models present) if (!m_Parameters.m_FiberModelList.empty()) { std::vector< double* > buffers; for (unsigned int i=0; iGetBufferPointer()); #pragma omp parallel for for( int i=0; iGetAbortGenerateData()) continue; float fiberWeight = m_FiberBundleTransformed->GetFiberWeight(i); int numPoints = -1; std::vector< itk::Vector > points_copy; #pragma omp critical { vtkCell* cell = fiberPolyData->GetCell(i); numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j))); } if (numPoints<2) continue; double seg_volume = fiberWeight*itk::Math::pi*m_mmRadius*m_mmRadius; for( int j=0; jGetAbortGenerateData()) { j=numPoints; continue; } itk::Vector v = points_copy.at(j); itk::Vector dir = points_copy.at(j+1)-v; if ( dir.GetSquaredNorm()<0.0001 || dir[0]!=dir[0] || dir[1]!=dir[1] || dir[2]!=dir[2] ) continue; dir.Normalize(); itk::Point startVertex = points_copy.at(j); itk::Index<3> startIndex; itk::ContinuousIndex startIndexCont; m_TransformedMaskImage->TransformPhysicalPointToIndex(startVertex, startIndex); m_TransformedMaskImage->TransformPhysicalPointToContinuousIndex(startVertex, startIndexCont); itk::Point endVertex = points_copy.at(j+1); itk::Index<3> endIndex; itk::ContinuousIndex endIndexCont; m_TransformedMaskImage->TransformPhysicalPointToIndex(endVertex, endIndex); m_TransformedMaskImage->TransformPhysicalPointToContinuousIndex(endVertex, endIndexCont); std::vector< std::pair< itk::Index<3>, double > > segments = mitk::imv::IntersectImage(m_WorkingSpacing, startIndex, endIndex, startIndexCont, endIndexCont); // generate signal for each fiber compartment for (int k=0; kSimulateMeasurement(g, dir)*seg_volume; for (std::pair< itk::Index<3>, double > seg : segments) { if (!m_TransformedMaskImage->GetLargestPossibleRegion().IsInside(seg.first) || m_TransformedMaskImage->GetPixel(seg.first)<=0) continue; double seg_signal = seg.second*signal_add; unsigned int linear_index = g + num_gradients*seg.first[0] + num_gradients*image_size_x*seg.first[1] + num_gradients*image_size_x*region_size_y*seg.first[2]; // update dMRI volume #pragma omp atomic buffers[k][linear_index] += seg_signal; // update fiber volume image if (k==0) { linear_index = seg.first[0] + image_size_x*seg.first[1] + image_size_x*region_size_y*seg.first[2]; #pragma omp atomic intraAxBuffer[linear_index] += seg.second*seg_volume; double vol = intraAxBuffer[linear_index]; if (vol>maxVolume) { maxVolume = vol; } } } } } #pragma omp critical { // progress report ++disp; unsigned long newTick = 50*disp.count()/disp.expected_count(); for (unsigned int tick = 0; tick<(newTick-lastTick); ++tick) PrintToLog("*", false, false, false); lastTick = newTick; } } } // axon radius not manually defined --> set fullest voxel (maxVolume) to full fiber voxel double density_correctiony_global = 1.0; if (m_Parameters.m_SignalGen.m_AxonRadius<0.0001) density_correctiony_global = m_VoxelVolume/maxVolume; // generate non-fiber signal ImageRegionIterator it3(m_TransformedMaskImage, m_TransformedMaskImage->GetLargestPossibleRegion()); while(!it3.IsAtEnd()) { if (it3.Get()>0) { DoubleDwiType::IndexType index = it3.GetIndex(); double iAxVolume = intraAxonalVolumeImage->GetPixel(index); // get non-transformed point (remove headmotion tranformation) // this point lives in the volume fraction image space itk::Point volume_fraction_point; if ( m_Parameters.m_SignalGen.m_DoAddMotion ) volume_fraction_point = GetMovedPoint(index, false); else m_TransformedMaskImage->TransformIndexToPhysicalPoint(index, volume_fraction_point); if (m_Parameters.m_SignalGen.m_DoDisablePartialVolume) { if (iAxVolume>0.0001) // scale fiber compartment to voxel { DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(index); pix[g] *= m_VoxelVolume/iAxVolume; m_CompartmentImages.at(0)->SetPixel(index, pix); if (g==0) m_VolumeFractions.at(0)->SetPixel(index, 1); } else { DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(index); pix[g] = 0; m_CompartmentImages.at(0)->SetPixel(index, pix); SimulateExtraAxonalSignal(index, volume_fraction_point, 0, g); } } else { // manually defined axon radius and voxel overflow --> rescale to voxel volume if ( m_Parameters.m_SignalGen.m_AxonRadius>=0.0001 && iAxVolume>m_VoxelVolume ) { for (int i=0; iGetPixel(index); pix[g] *= m_VoxelVolume/iAxVolume; m_CompartmentImages.at(i)->SetPixel(index, pix); } iAxVolume = m_VoxelVolume; } // if volume fraction image is set use it, otherwise use global scaling factor double density_correction_voxel = density_correctiony_global; if ( m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()!=nullptr && iAxVolume>0.0001 ) { m_DoubleInterpolator->SetInputImage(m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()); double volume_fraction = mitk::imv::GetImageValue(volume_fraction_point, true, m_DoubleInterpolator); if (volume_fraction<0) mitkThrow() << "Volume fraction image (index 1) contains negative values (intra-axonal compartment)!"; density_correction_voxel = m_VoxelVolume*volume_fraction/iAxVolume; // remove iAxVolume sclaing and scale to volume_fraction } else if (m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()!=nullptr) density_correction_voxel = 0.0; // adjust intra-axonal compartment volume by density correction factor DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(index); pix[g] *= density_correction_voxel; m_CompartmentImages.at(0)->SetPixel(index, pix); // normalize remaining fiber volume fractions (they are rescaled in SimulateExtraAxonalSignal) if (iAxVolume>0.0001) { for (int i=1; iGetPixel(index); pix[g] /= iAxVolume; m_CompartmentImages.at(i)->SetPixel(index, pix); } } else { for (int i=1; iGetPixel(index); pix[g] = 0; m_CompartmentImages.at(i)->SetPixel(index, pix); } } iAxVolume = density_correction_voxel*iAxVolume; // new intra-axonal volume = old intra-axonal volume * correction factor // simulate other compartments SimulateExtraAxonalSignal(index, volume_fraction_point, iAxVolume, g); } } ++it3; } } PrintToLog("\n", false); } if (this->GetAbortGenerateData()) { PrintToLog("\n", false, false); PrintToLog("Simulation aborted"); return; } DoubleDwiType::Pointer doubleOutImage; double signalScale = m_Parameters.m_SignalGen.m_SignalScale; if ( m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition ) // do k-space stuff { PrintToLog("\n", false, false); PrintToLog("Simulating k-space acquisition using " +boost::lexical_cast(m_Parameters.m_SignalGen.m_NumberOfCoils) +" coil(s)"); switch (m_Parameters.m_SignalGen.m_AcquisitionType) { case SignalGenerationParameters::SingleShotEpi: { PrintToLog("Acquisition type: single shot EPI", false); break; } case SignalGenerationParameters::SpinEcho: { PrintToLog("Acquisition type: classic spin echo with cartesian k-space trajectory", false); break; } default: { PrintToLog("Acquisition type: single shot EPI", false); break; } } if (m_Parameters.m_SignalGen.m_DoSimulateRelaxation) PrintToLog("Simulating signal relaxation", false); if (m_Parameters.m_SignalGen.m_NoiseVariance>0 && m_Parameters.m_Misc.m_DoAddNoise) PrintToLog("Simulating complex Gaussian noise: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_NoiseVariance), false); if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull() && m_Parameters.m_Misc.m_DoAddDistortions) PrintToLog("Simulating distortions", false); if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging) PrintToLog("Simulating ringing artifacts", false); if (m_Parameters.m_Misc.m_DoAddEddyCurrents && m_Parameters.m_SignalGen.m_EddyStrength>0) PrintToLog("Simulating eddy currents: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_EddyStrength), false); if (m_Parameters.m_Misc.m_DoAddSpikes && m_Parameters.m_SignalGen.m_Spikes>0) PrintToLog("Simulating spikes: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_Spikes), false); if (m_Parameters.m_Misc.m_DoAddAliasing && m_Parameters.m_SignalGen.m_CroppingFactor<1.0) PrintToLog("Simulating aliasing: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_CroppingFactor), false); if (m_Parameters.m_Misc.m_DoAddGhosts && m_Parameters.m_SignalGen.m_KspaceLineOffset>0) PrintToLog("Simulating ghosts: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_KspaceLineOffset), false); doubleOutImage = SimulateKspaceAcquisition(m_CompartmentImages); signalScale = 1; // already scaled in SimulateKspaceAcquisition() } else // don't do k-space stuff, just sum compartments { PrintToLog("Summing compartments"); doubleOutImage = m_CompartmentImages.at(0); for (unsigned int i=1; i::New(); adder->SetInput1(doubleOutImage); adder->SetInput2(m_CompartmentImages.at(i)); adder->Update(); doubleOutImage = adder->GetOutput(); } } if (this->GetAbortGenerateData()) { PrintToLog("\n", false, false); PrintToLog("Simulation aborted"); return; } PrintToLog("Finalizing image"); if (m_Parameters.m_SignalGen.m_DoAddDrift && m_Parameters.m_SignalGen.m_Drift>0.0) PrintToLog("Adding signal drift: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_Drift), false); if (signalScale>1) PrintToLog("Scaling signal", false); if (m_Parameters.m_NoiseModel) PrintToLog("Adding noise: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_NoiseVariance), false); unsigned int window = 0; unsigned int min = itk::NumericTraits::max(); ImageRegionIterator it4 (m_OutputImage, m_OutputImage->GetLargestPossibleRegion()); DoubleDwiType::PixelType signal; signal.SetSize(m_Parameters.m_SignalGen.GetNumVolumes()); boost::progress_display disp2(m_OutputImage->GetLargestPossibleRegion().GetNumberOfPixels()); PrintToLog("0% 10 20 30 40 50 60 70 80 90 100%", false, true, false); PrintToLog("|----|----|----|----|----|----|----|----|----|----|\n*", false, false, false); int lastTick = 0; while(!it4.IsAtEnd()) { if (this->GetAbortGenerateData()) { PrintToLog("\n", false, false); PrintToLog("Simulation aborted"); return; } ++disp2; unsigned long newTick = 50*disp2.count()/disp2.expected_count(); for (unsigned long tick = 0; tick<(newTick-lastTick); tick++) PrintToLog("*", false, false, false); lastTick = newTick; typename OutputImageType::IndexType index = it4.GetIndex(); signal = doubleOutImage->GetPixel(index)*signalScale; for (unsigned int i=0; iAddNoise(signal); for (unsigned int i=0; i0) signal[i] = floor(signal[i]+0.5); else signal[i] = ceil(signal[i]-0.5); if ( (!m_Parameters.m_SignalGen.IsBaselineIndex(i) || signal.Size()==1) && signal[i]>window) window = signal[i]; if ( (!m_Parameters.m_SignalGen.IsBaselineIndex(i) || signal.Size()==1) && signal[i]SetNthOutput(0, m_OutputImage); PrintToLog("\n", false); PrintToLog("Finished simulation"); m_TimeProbe.Stop(); if (m_Parameters.m_SignalGen.m_DoAddMotion) { PrintToLog("\nHead motion log:", false); PrintToLog(m_MotionLog, false, false); } if (m_Parameters.m_Misc.m_DoAddSpikes && m_Parameters.m_SignalGen.m_Spikes>0) { PrintToLog("\nSpike log:", false); PrintToLog(m_SpikeLog, false, false); } if (m_Logfile.is_open()) m_Logfile.close(); } template< class PixelType > void TractsToDWIImageFilter< PixelType >::PrintToLog(std::string m, bool addTime, bool linebreak, bool stdOut) { // timestamp if (addTime) { m_Logfile << this->GetTime() << " > "; m_StatusText += this->GetTime() + " > "; if (stdOut) std::cout << this->GetTime() << " > "; } // message if (m_Logfile.is_open()) m_Logfile << m; m_StatusText += m; if (stdOut) std::cout << m; // new line if (linebreak) { if (m_Logfile.is_open()) m_Logfile << "\n"; m_StatusText += "\n"; if (stdOut) std::cout << "\n"; } m_Logfile.flush(); } template< class PixelType > void TractsToDWIImageFilter< PixelType >::SimulateMotion(int g) { if ( m_Parameters.m_SignalGen.m_DoAddMotion && m_Parameters.m_SignalGen.m_DoRandomizeMotion && g>0 && m_Parameters.m_SignalGen.m_MotionVolumes[g-1]) { // The last volume was randomly moved, so we have to reset to fiberbundle and the mask. // Without motion or with linear motion, we keep the last position --> no reset. m_FiberBundleTransformed = m_FiberBundle->GetDeepCopy(); if (m_MaskImageSet) { auto duplicator = itk::ImageDuplicator::New(); duplicator->SetInputImage(m_Parameters.m_SignalGen.m_MaskImage); duplicator->Update(); m_TransformedMaskImage = duplicator->GetOutput(); } } VectorType rotation; VectorType translation; // is motion artifact enabled? // is the current volume g affected by motion? if ( m_Parameters.m_SignalGen.m_DoAddMotion && m_Parameters.m_SignalGen.m_MotionVolumes[g] && g(m_Parameters.m_SignalGen.GetNumVolumes()) ) { // adjust motion transforms if ( m_Parameters.m_SignalGen.m_DoRandomizeMotion ) { // randomly rotation[0] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[0]*2) -m_Parameters.m_SignalGen.m_Rotation[0]; rotation[1] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[1]*2) -m_Parameters.m_SignalGen.m_Rotation[1]; rotation[2] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[2]*2) -m_Parameters.m_SignalGen.m_Rotation[2]; translation[0] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[0]*2) -m_Parameters.m_SignalGen.m_Translation[0]; translation[1] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[1]*2) -m_Parameters.m_SignalGen.m_Translation[1]; translation[2] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[2]*2) -m_Parameters.m_SignalGen.m_Translation[2]; m_FiberBundleTransformed->TransformFibers(rotation[0], rotation[1], rotation[2], translation[0], translation[1], translation[2]); } else { // linearly rotation = m_Parameters.m_SignalGen.m_Rotation / m_NumMotionVolumes; translation = m_Parameters.m_SignalGen.m_Translation / m_NumMotionVolumes; m_MotionCounter++; m_FiberBundleTransformed->TransformFibers(rotation[0], rotation[1], rotation[2], translation[0], translation[1], translation[2]); rotation *= m_MotionCounter; translation *= m_MotionCounter; } MatrixType rotationMatrix = mitk::imv::GetRotationMatrixItk(rotation[0], rotation[1], rotation[2]); MatrixType rotationMatrixInv = mitk::imv::GetRotationMatrixItk(-rotation[0], -rotation[1], -rotation[2]); m_Rotations.push_back(rotationMatrix); m_RotationsInv.push_back(rotationMatrixInv); m_Translations.push_back(translation); // move mask image accoring to new transform if (m_MaskImageSet) { ImageRegionIterator maskIt(m_UpsampledMaskImage, m_UpsampledMaskImage->GetLargestPossibleRegion()); m_TransformedMaskImage->FillBuffer(0); while(!maskIt.IsAtEnd()) { if (maskIt.Get()<=0) { ++maskIt; continue; } DoubleDwiType::IndexType index = maskIt.GetIndex(); m_TransformedMaskImage->TransformPhysicalPointToIndex(GetMovedPoint(index, true), index); if (m_TransformedMaskImage->GetLargestPossibleRegion().IsInside(index)) m_TransformedMaskImage->SetPixel(index, 100); ++maskIt; } } } else { if (m_Parameters.m_SignalGen.m_DoAddMotion && !m_Parameters.m_SignalGen.m_DoRandomizeMotion && g>0) { rotation = m_Parameters.m_SignalGen.m_Rotation / m_NumMotionVolumes; rotation *= m_MotionCounter; m_Rotations.push_back(m_Rotations.back()); m_RotationsInv.push_back(m_RotationsInv.back()); m_Translations.push_back(m_Translations.back()); } else { rotation.Fill(0.0); VectorType translation; translation.Fill(0.0); MatrixType rotation_matrix; rotation_matrix.SetIdentity(); m_Rotations.push_back(rotation_matrix); m_RotationsInv.push_back(rotation_matrix); m_Translations.push_back(translation); } } if (m_Parameters.m_SignalGen.m_DoAddMotion) { m_MotionLog += boost::lexical_cast(g) + " rotation: " + boost::lexical_cast(rotation[0]) + "," + boost::lexical_cast(rotation[1]) + "," + boost::lexical_cast(rotation[2]) + ";"; m_MotionLog += " translation: " + boost::lexical_cast(m_Translations.back()[0]) + "," + boost::lexical_cast(m_Translations.back()[1]) + "," + boost::lexical_cast(m_Translations.back()[2]) + "\n"; } } template< class PixelType > itk::Point TractsToDWIImageFilter< PixelType >::GetMovedPoint(itk::Index<3>& index, bool forward) { itk::Point transformed_point; float tx = m_Translations.back()[0]; float ty = m_Translations.back()[1]; float tz = m_Translations.back()[2]; if (forward) { m_UpsampledMaskImage->TransformIndexToPhysicalPoint(index, transformed_point); m_FiberBundle->TransformPoint<>(transformed_point, m_Rotations.back(), tx, ty, tz); } else { tx *= -1; ty *= -1; tz *= -1; m_TransformedMaskImage->TransformIndexToPhysicalPoint(index, transformed_point); m_FiberBundle->TransformPoint<>(transformed_point, m_RotationsInv.back(), tx, ty, tz); } return transformed_point; } template< class PixelType > void TractsToDWIImageFilter< PixelType >:: SimulateExtraAxonalSignal(ItkUcharImgType::IndexType& index, itk::Point& volume_fraction_point, double intraAxonalVolume, int g) { int numFiberCompartments = m_Parameters.m_FiberModelList.size(); int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size(); if (m_Parameters.m_SignalGen.m_DoDisablePartialVolume) { // simulate signal for largest non-fiber compartment int max_compartment_index = 0; double max_fraction = 0; if (numNonFiberCompartments>1) { for (int i=0; iSetInputImage(m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()); double compartment_fraction = mitk::imv::GetImageValue(volume_fraction_point, true, m_DoubleInterpolator); if (compartment_fraction<0) mitkThrow() << "Volume fraction image (index " << i << ") contains values less than zero!"; if (compartment_fraction>max_fraction) { max_fraction = compartment_fraction; max_compartment_index = i; } } } DoubleDwiType::Pointer doubleDwi = m_CompartmentImages.at(max_compartment_index+numFiberCompartments); DoubleDwiType::PixelType pix = doubleDwi->GetPixel(index); pix[g] += m_Parameters.m_NonFiberModelList[max_compartment_index]->SimulateMeasurement(g, m_NullDir)*m_VoxelVolume; doubleDwi->SetPixel(index, pix); if (g==0) m_VolumeFractions.at(max_compartment_index+numFiberCompartments)->SetPixel(index, 1); } else { std::vector< double > fractions; if (g==0) m_VolumeFractions.at(0)->SetPixel(index, intraAxonalVolume/m_VoxelVolume); double extraAxonalVolume = m_VoxelVolume-intraAxonalVolume; // non-fiber volume if (extraAxonalVolume<0) { if (extraAxonalVolume<-0.001) MITK_ERROR << "Corrupted intra-axonal signal voxel detected. Fiber volume larger voxel volume! " << m_VoxelVolume << "<" << intraAxonalVolume; extraAxonalVolume = 0; } double interAxonalVolume = 0; if (numFiberCompartments>1) interAxonalVolume = extraAxonalVolume * intraAxonalVolume/m_VoxelVolume; // inter-axonal fraction of non fiber compartment double nonFiberVolume = extraAxonalVolume - interAxonalVolume; // rest of compartment if (nonFiberVolume<0) { if (nonFiberVolume<-0.001) MITK_ERROR << "Corrupted signal voxel detected. Fiber volume larger voxel volume!"; nonFiberVolume = 0; interAxonalVolume = extraAxonalVolume; } double compartmentSum = intraAxonalVolume; fractions.push_back(intraAxonalVolume/m_VoxelVolume); // rescale extra-axonal fiber signal for (int i=1; iGetVolumeFractionImage()!=nullptr) { m_DoubleInterpolator->SetInputImage(m_Parameters.m_FiberModelList[i]->GetVolumeFractionImage()); interAxonalVolume = mitk::imv::GetImageValue(volume_fraction_point, true, m_DoubleInterpolator)*m_VoxelVolume; if (interAxonalVolume<0) mitkThrow() << "Volume fraction image (index " << i+1 << ") contains negative values!"; } DoubleDwiType::PixelType pix = m_CompartmentImages.at(i)->GetPixel(index); pix[g] *= interAxonalVolume; m_CompartmentImages.at(i)->SetPixel(index, pix); compartmentSum += interAxonalVolume; fractions.push_back(interAxonalVolume/m_VoxelVolume); if (g==0) m_VolumeFractions.at(i)->SetPixel(index, interAxonalVolume/m_VoxelVolume); } for (int i=0; iGetVolumeFractionImage()!=nullptr) { m_DoubleInterpolator->SetInputImage(m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()); volume = mitk::imv::GetImageValue(volume_fraction_point, true, m_DoubleInterpolator)*m_VoxelVolume; if (volume<0) mitkThrow() << "Volume fraction image (index " << numFiberCompartments+i+1 << ") contains negative values (non-fiber compartment)!"; if (m_UseRelativeNonFiberVolumeFractions) volume *= nonFiberVolume/m_VoxelVolume; } DoubleDwiType::PixelType pix = m_CompartmentImages.at(i+numFiberCompartments)->GetPixel(index); pix[g] += m_Parameters.m_NonFiberModelList[i]->SimulateMeasurement(g, m_NullDir)*volume; m_CompartmentImages.at(i+numFiberCompartments)->SetPixel(index, pix); compartmentSum += volume; fractions.push_back(volume/m_VoxelVolume); if (g==0) m_VolumeFractions.at(i+numFiberCompartments)->SetPixel(index, volume/m_VoxelVolume); } if (compartmentSum/m_VoxelVolume>1.05) { MITK_ERROR << "Compartments do not sum to 1 in voxel " << index << " (" << compartmentSum/m_VoxelVolume << ")"; for (auto val : fractions) MITK_ERROR << val; } } } template< class PixelType > itk::Vector TractsToDWIImageFilter< PixelType >::GetItkVector(double point[3]) { itk::Vector itkVector; itkVector[0] = point[0]; itkVector[1] = point[1]; itkVector[2] = point[2]; return itkVector; } template< class PixelType > vnl_vector_fixed TractsToDWIImageFilter< PixelType >::GetVnlVector(double point[3]) { vnl_vector_fixed vnlVector; vnlVector[0] = point[0]; vnlVector[1] = point[1]; vnlVector[2] = point[2]; return vnlVector; } template< class PixelType > vnl_vector_fixed TractsToDWIImageFilter< PixelType >::GetVnlVector(Vector& vector) { vnl_vector_fixed vnlVector; vnlVector[0] = vector[0]; vnlVector[1] = vector[1]; vnlVector[2] = vector[2]; return vnlVector; } template< class PixelType > double TractsToDWIImageFilter< PixelType >::RoundToNearest(double num) { return (num > 0.0) ? floor(num + 0.5) : ceil(num - 0.5); } template< class PixelType > std::string TractsToDWIImageFilter< PixelType >::GetTime() { m_TimeProbe.Stop(); unsigned long total = RoundToNearest(m_TimeProbe.GetTotal()); unsigned long hours = total/3600; unsigned long minutes = (total%3600)/60; unsigned long seconds = total%60; std::string out = ""; out.append(boost::lexical_cast(hours)); out.append(":"); out.append(boost::lexical_cast(minutes)); out.append(":"); out.append(boost::lexical_cast(seconds)); m_TimeProbe.Start(); return out; } } diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.cpp index a0d157612d..81d6933277 100644 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.cpp +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.cpp @@ -1,1068 +1,1071 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #define RAPIDXML_NO_EXCEPTIONS #include #include #include #include #include #include #include #include #include mitk::FiberfoxParameters::FiberfoxParameters() : m_NoiseModel(nullptr) { mitk::StickModel* model_aniso = new mitk::StickModel(); model_aniso->m_CompartmentId = 1; m_FiberModelList.push_back(model_aniso); mitk::BallModel* model_iso = new mitk::BallModel(); model_iso->m_CompartmentId = 3; m_NonFiberModelList.push_back(model_iso); } mitk::FiberfoxParameters::FiberfoxParameters(const mitk::FiberfoxParameters& params) : m_NoiseModel(nullptr) { m_FiberGen = params.m_FiberGen; m_SignalGen = params.m_SignalGen; m_Misc = params.m_Misc; if (params.m_NoiseModel!=nullptr) { if (dynamic_cast*>(params.m_NoiseModel.get())) m_NoiseModel = std::make_shared< mitk::RicianNoiseModel<> >(); else if (dynamic_cast*>(params.m_NoiseModel.get())) m_NoiseModel = std::make_shared< mitk::ChiSquareNoiseModel<> >(); m_NoiseModel->SetNoiseVariance(params.m_NoiseModel->GetNoiseVariance()); } for (unsigned int i=0; i* outModel = nullptr; mitk::DiffusionSignalModel<>* signalModel = nullptr; if (i*>(signalModel)) outModel = new mitk::StickModel<>(dynamic_cast*>(signalModel)); else if (dynamic_cast*>(signalModel)) outModel = new mitk::TensorModel<>(dynamic_cast*>(signalModel)); else if (dynamic_cast*>(signalModel)) outModel = new mitk::RawShModel<>(dynamic_cast*>(signalModel)); else if (dynamic_cast*>(signalModel)) outModel = new mitk::BallModel<>(dynamic_cast*>(signalModel)); else if (dynamic_cast*>(signalModel)) outModel = new mitk::AstroStickModel<>(dynamic_cast*>(signalModel)); else if (dynamic_cast*>(signalModel)) outModel = new mitk::DotModel<>(dynamic_cast*>(signalModel)); if (i theta; theta.set_size(NPoints); vnl_vector phi; phi.set_size(NPoints); double C = sqrt(4*itk::Math::pi); phi(0) = 0.0; phi(NPoints-1) = 0.0; for(int i=0; i0 && i mitk::SignalGenerationParameters::GetBaselineIndices() { std::vector< int > result; for( unsigned int i=0; im_GradientDirections.size(); i++) if (m_GradientDirections.at(i).GetNorm()<0.0001) result.push_back(i); return result; } unsigned int mitk::SignalGenerationParameters::GetFirstBaselineIndex() { for( unsigned int i=0; im_GradientDirections.size(); i++) if (m_GradientDirections.at(i).GetNorm()<0.0001) return i; return -1; } bool mitk::SignalGenerationParameters::IsBaselineIndex(unsigned int idx) { if (m_GradientDirections.size()>idx && m_GradientDirections.at(idx).GetNorm()<0.0001) return true; return false; } unsigned int mitk::SignalGenerationParameters::GetNumWeightedVolumes() { return m_NumGradients; } unsigned int mitk::SignalGenerationParameters::GetNumBaselineVolumes() { return m_NumBaseline; } unsigned int mitk::SignalGenerationParameters::GetNumVolumes() { return m_GradientDirections.size(); } mitk::SignalGenerationParameters::GradientListType mitk::SignalGenerationParameters::GetGradientDirections() { return m_GradientDirections; } mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer mitk::SignalGenerationParameters::GetItkGradientContainer() { int c = 0; mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer out = mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::New(); for (auto g : m_GradientDirections) { mitk::DiffusionPropertyHelper::GradientDirectionType vnl_dir; vnl_dir[0] = g[0]; vnl_dir[1] = g[1]; vnl_dir[2] = g[2]; out->InsertElement(c, vnl_dir); ++c; } return out; } mitk::SignalGenerationParameters::GradientType mitk::SignalGenerationParameters::GetGradientDirection(unsigned int i) { return m_GradientDirections.at(i); } void mitk::SignalGenerationParameters::SetNumWeightedVolumes(int numGradients) { m_NumGradients = numGradients; GenerateGradientHalfShell(); } std::vector< int > mitk::SignalGenerationParameters::GetBvalues() { std::vector< int > bVals; for( GradientType g : m_GradientDirections) { float norm = g.GetNorm(); int bVal = std::round(norm*norm*m_Bvalue); if ( std::find(bVals.begin(), bVals.end(), bVal) == bVals.end() ) bVals.push_back(bVal); } return bVals; } double mitk::SignalGenerationParameters::GetBvalue() { return m_Bvalue; } void mitk::SignalGenerationParameters::SetGradienDirections(GradientListType gradientList) { m_GradientDirections = gradientList; m_NumGradients = 0; m_NumBaseline = 0; for( unsigned int i=0; im_GradientDirections.size(); i++) { float norm = m_GradientDirections.at(i).GetNorm(); if (norm>0.0001) m_NumGradients++; else m_NumBaseline++; } } void mitk::SignalGenerationParameters::SetGradienDirections(mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradientList) { m_NumGradients = 0; m_NumBaseline = 0; m_GradientDirections.clear(); for( unsigned int i=0; iSize(); i++) { GradientType g; g[0] = gradientList->at(i)[0]; g[1] = gradientList->at(i)[1]; g[2] = gradientList->at(i)[2]; m_GradientDirections.push_back(g); float norm = m_GradientDirections.at(i).GetNorm(); if (norm>0.0001) m_NumGradients++; else m_NumBaseline++; } } void mitk::SignalGenerationParameters::ApplyDirectionMatrix() { auto imageRotationMatrix = m_ImageDirection.GetVnlMatrix(); GradientListType rotated_gradients; for(auto g : m_GradientDirections) { vnl_vector vec = g.GetVnlVector(); vec = vec.pre_multiply(imageRotationMatrix); GradientType g2; g2[0] = vec[0]; g2[1] = vec[1]; g2[2] = vec[2]; rotated_gradients.push_back(g2); } m_GradientDirections = rotated_gradients; } void mitk::FiberfoxParameters::ApplyDirectionMatrix() { m_SignalGen.ApplyDirectionMatrix(); UpdateSignalModels(); } void mitk::FiberfoxParameters::SaveParameters(std::string filename) { if(filename.empty()) return; if(".ffp"!=filename.substr(filename.size()-4, 4)) filename += ".ffp"; const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, nullptr ); if ( locale.compare(currLocale)!=0 ) { try { setlocale(LC_ALL, locale.c_str()); } catch(...) { MITK_INFO << "Could not set locale " << locale; } } boost::property_tree::ptree parameters; // fiber generation parameters parameters.put("fiberfox.fibers.distribution", m_FiberGen.m_Distribution); parameters.put("fiberfox.fibers.variance", m_FiberGen.m_Variance); parameters.put("fiberfox.fibers.density", m_FiberGen.m_Density); parameters.put("fiberfox.fibers.spline.sampling", m_FiberGen.m_Sampling); parameters.put("fiberfox.fibers.spline.tension", m_FiberGen.m_Tension); parameters.put("fiberfox.fibers.spline.continuity", m_FiberGen.m_Continuity); parameters.put("fiberfox.fibers.spline.bias", m_FiberGen.m_Bias); parameters.put("fiberfox.fibers.rotation.x", m_FiberGen.m_Rotation[0]); parameters.put("fiberfox.fibers.rotation.y", m_FiberGen.m_Rotation[1]); parameters.put("fiberfox.fibers.rotation.z", m_FiberGen.m_Rotation[2]); parameters.put("fiberfox.fibers.translation.x", m_FiberGen.m_Translation[0]); parameters.put("fiberfox.fibers.translation.y", m_FiberGen.m_Translation[1]); parameters.put("fiberfox.fibers.translation.z", m_FiberGen.m_Translation[2]); parameters.put("fiberfox.fibers.scale.x", m_FiberGen.m_Scale[0]); parameters.put("fiberfox.fibers.scale.y", m_FiberGen.m_Scale[1]); parameters.put("fiberfox.fibers.scale.z", m_FiberGen.m_Scale[2]); // image generation parameters parameters.put("fiberfox.image.basic.size.x", m_SignalGen.m_ImageRegion.GetSize(0)); parameters.put("fiberfox.image.basic.size.y", m_SignalGen.m_ImageRegion.GetSize(1)); parameters.put("fiberfox.image.basic.size.z", m_SignalGen.m_ImageRegion.GetSize(2)); parameters.put("fiberfox.image.basic.spacing.x", m_SignalGen.m_ImageSpacing[0]); parameters.put("fiberfox.image.basic.spacing.y", m_SignalGen.m_ImageSpacing[1]); parameters.put("fiberfox.image.basic.spacing.z", m_SignalGen.m_ImageSpacing[2]); parameters.put("fiberfox.image.basic.origin.x", m_SignalGen.m_ImageOrigin[0]); parameters.put("fiberfox.image.basic.origin.y", m_SignalGen.m_ImageOrigin[1]); parameters.put("fiberfox.image.basic.origin.z", m_SignalGen.m_ImageOrigin[2]); parameters.put("fiberfox.image.basic.direction.d1", m_SignalGen.m_ImageDirection[0][0]); parameters.put("fiberfox.image.basic.direction.d2", m_SignalGen.m_ImageDirection[0][1]); parameters.put("fiberfox.image.basic.direction.d3", m_SignalGen.m_ImageDirection[0][2]); parameters.put("fiberfox.image.basic.direction.d4", m_SignalGen.m_ImageDirection[1][0]); parameters.put("fiberfox.image.basic.direction.d5", m_SignalGen.m_ImageDirection[1][1]); parameters.put("fiberfox.image.basic.direction.d6", m_SignalGen.m_ImageDirection[1][2]); parameters.put("fiberfox.image.basic.direction.d7", m_SignalGen.m_ImageDirection[2][0]); parameters.put("fiberfox.image.basic.direction.d8", m_SignalGen.m_ImageDirection[2][1]); parameters.put("fiberfox.image.basic.direction.d9", m_SignalGen.m_ImageDirection[2][2]); mitk::gradients::WriteBvalsBvecs(filename+".bvals", filename+".bvecs", m_SignalGen.GetItkGradientContainer(), m_SignalGen.m_Bvalue); parameters.put("fiberfox.image.acquisitiontype", m_SignalGen.m_AcquisitionType); parameters.put("fiberfox.image.coilsensitivityprofile", m_SignalGen.m_CoilSensitivityProfile); parameters.put("fiberfox.image.numberofcoils", m_SignalGen.m_NumberOfCoils); parameters.put("fiberfox.image.reversephase", m_SignalGen.m_ReversePhase); parameters.put("fiberfox.image.partialfourier", m_SignalGen.m_PartialFourier); parameters.put("fiberfox.image.noisevariance", m_SignalGen.m_NoiseVariance); parameters.put("fiberfox.image.trep", m_SignalGen.m_tRep); parameters.put("fiberfox.image.signalScale", m_SignalGen.m_SignalScale); parameters.put("fiberfox.image.tEcho", m_SignalGen.m_tEcho); parameters.put("fiberfox.image.tLine", m_SignalGen.m_tLine); parameters.put("fiberfox.image.tInhom", m_SignalGen.m_tInhom); parameters.put("fiberfox.image.simulatekspace", m_SignalGen.m_SimulateKspaceAcquisition); parameters.put("fiberfox.image.axonRadius", m_SignalGen.m_AxonRadius); parameters.put("fiberfox.image.doSimulateRelaxation", m_SignalGen.m_DoSimulateRelaxation); parameters.put("fiberfox.image.doDisablePartialVolume", m_SignalGen.m_DoDisablePartialVolume); parameters.put("fiberfox.image.artifacts.spikesnum", m_SignalGen.m_Spikes); parameters.put("fiberfox.image.artifacts.spikesscale", m_SignalGen.m_SpikeAmplitude); parameters.put("fiberfox.image.artifacts.kspaceLineOffset", m_SignalGen.m_KspaceLineOffset); parameters.put("fiberfox.image.artifacts.eddyStrength", m_SignalGen.m_EddyStrength); parameters.put("fiberfox.image.artifacts.eddyTau", m_SignalGen.m_Tau); parameters.put("fiberfox.image.artifacts.aliasingfactor", m_SignalGen.m_CroppingFactor); parameters.put("fiberfox.image.artifacts.drift", m_SignalGen.m_Drift); parameters.put("fiberfox.image.artifacts.doAddMotion", m_SignalGen.m_DoAddMotion); parameters.put("fiberfox.image.artifacts.randomMotion", m_SignalGen.m_DoRandomizeMotion); parameters.put("fiberfox.image.artifacts.translation0", m_SignalGen.m_Translation[0]); parameters.put("fiberfox.image.artifacts.translation1", m_SignalGen.m_Translation[1]); parameters.put("fiberfox.image.artifacts.translation2", m_SignalGen.m_Translation[2]); parameters.put("fiberfox.image.artifacts.rotation0", m_SignalGen.m_Rotation[0]); parameters.put("fiberfox.image.artifacts.rotation1", m_SignalGen.m_Rotation[1]); parameters.put("fiberfox.image.artifacts.rotation2", m_SignalGen.m_Rotation[2]); parameters.put("fiberfox.image.artifacts.motionvolumes", m_Misc.m_MotionVolumesBox); parameters.put("fiberfox.image.artifacts.addringing", m_SignalGen.m_DoAddGibbsRinging); parameters.put("fiberfox.image.artifacts.addnoise", m_Misc.m_DoAddNoise); parameters.put("fiberfox.image.artifacts.addghosts", m_Misc.m_DoAddGhosts); parameters.put("fiberfox.image.artifacts.addaliasing", m_Misc.m_DoAddAliasing); parameters.put("fiberfox.image.artifacts.addspikes", m_Misc.m_DoAddSpikes); parameters.put("fiberfox.image.artifacts.addeddycurrents", m_Misc.m_DoAddEddyCurrents); parameters.put("fiberfox.image.artifacts.doAddDistortions", m_Misc.m_DoAddDistortions); parameters.put("fiberfox.image.artifacts.doAddDrift", m_SignalGen.m_DoAddDrift); parameters.put("fiberfox.image.outputvolumefractions", m_Misc.m_CheckOutputVolumeFractionsBox); parameters.put("fiberfox.image.showadvanced", m_Misc.m_CheckAdvancedSignalOptionsBox); parameters.put("fiberfox.image.signalmodelstring", m_Misc.m_SignalModelString); parameters.put("fiberfox.image.artifactmodelstring", m_Misc.m_ArtifactModelString); parameters.put("fiberfox.image.outpath", m_Misc.m_OutputPath); parameters.put("fiberfox.fibers.realtime", m_Misc.m_CheckRealTimeFibersBox); parameters.put("fiberfox.fibers.showadvanced", m_Misc.m_CheckAdvancedFiberOptionsBox); parameters.put("fiberfox.fibers.constantradius", m_Misc.m_CheckConstantRadiusBox); parameters.put("fiberfox.fibers.includeFiducials", m_Misc.m_CheckIncludeFiducialsBox); if (m_NoiseModel!=nullptr) { parameters.put("fiberfox.image.artifacts.noisevariance", m_NoiseModel->GetNoiseVariance()); if (dynamic_cast*>(m_NoiseModel.get())) parameters.put("fiberfox.image.artifacts.noisetype", "rice"); else if (dynamic_cast*>(m_NoiseModel.get())) parameters.put("fiberfox.image.artifacts.noisetype", "chisquare"); } for (std::size_t i=0; i* signalModel = nullptr; if (i(i)+".type", "fiber"); } else { signalModel = m_NonFiberModelList.at(i-m_FiberModelList.size()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".type", "non-fiber"); } if (dynamic_cast*>(signalModel)) { mitk::StickModel<>* model = dynamic_cast*>(signalModel); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".model", "stick"); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".d", model->GetDiffusivity()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".t2", model->GetT2()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".t1", model->GetT1()); } else if (dynamic_cast*>(signalModel)) { mitk::TensorModel<>* model = dynamic_cast*>(signalModel); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".model", "tensor"); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".d1", model->GetDiffusivity1()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".d2", model->GetDiffusivity2()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".d3", model->GetDiffusivity3()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".t2", model->GetT2()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".t1", model->GetT1()); } else if (dynamic_cast*>(signalModel)) { mitk::RawShModel<>* model = dynamic_cast*>(signalModel); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".model", "prototype"); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".minFA", model->GetFaRange().first); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".maxFA", model->GetFaRange().second); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".minADC", model->GetAdcRange().first); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".maxADC", model->GetAdcRange().second); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".maxNumSamples", model->GetMaxNumKernels()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".numSamples", model->GetNumberOfKernels()); int shOrder = model->GetShOrder(); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".numCoeffs", (shOrder*shOrder + shOrder + 2)/2 + shOrder); for (unsigned int j=0; jGetNumberOfKernels(); j++) { vnl_vector< double > coeffs = model->GetCoefficients(j); for (unsigned int k=0; k(i)+".kernels."+boost::lexical_cast(j)+".coeffs."+boost::lexical_cast(k), coeffs[k]); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".kernels."+boost::lexical_cast(j)+".B0", model->GetBaselineSignal(j)); } } else if (dynamic_cast*>(signalModel)) { mitk::BallModel<>* model = dynamic_cast*>(signalModel); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".model", "ball"); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".d", model->GetDiffusivity()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".t2", model->GetT2()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".t1", model->GetT1()); } else if (dynamic_cast*>(signalModel)) { mitk::AstroStickModel<>* model = dynamic_cast*>(signalModel); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".model", "astrosticks"); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".d", model->GetDiffusivity()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".t2", model->GetT2()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".t1", model->GetT1()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".randomize", model->GetRandomizeSticks()); } else if (dynamic_cast*>(signalModel)) { mitk::DotModel<>* model = dynamic_cast*>(signalModel); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".model", "dot"); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".t2", model->GetT2()); parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".t1", model->GetT1()); } if (signalModel!=nullptr) { parameters.put("fiberfox.image.compartments.c"+boost::lexical_cast(i)+".ID", signalModel->m_CompartmentId); if (signalModel->GetVolumeFractionImage().IsNotNull()) { try{ itk::ImageFileWriter::Pointer writer = itk::ImageFileWriter::New(); writer->SetFileName(filename+"_VOLUME"+boost::lexical_cast(signalModel->m_CompartmentId)+".nii.gz"); writer->SetInput(signalModel->GetVolumeFractionImage()); writer->Update(); MITK_INFO << "Volume fraction image for compartment "+boost::lexical_cast(signalModel->m_CompartmentId)+" saved."; } catch(...) { } } } } boost::property_tree::xml_writer_settings writerSettings(' ', 2); boost::property_tree::xml_parser::write_xml(filename, parameters, std::locale(), writerSettings); try{ itk::ImageFileWriter::Pointer writer = itk::ImageFileWriter::New(); writer->SetFileName(filename+"_FMAP.nii.gz"); writer->SetInput(m_SignalGen.m_FrequencyMap); writer->Update(); } catch(...) { MITK_INFO << "No frequency map saved."; } try{ itk::ImageFileWriter::Pointer writer = itk::ImageFileWriter::New(); writer->SetFileName(filename+"_MASK.nii.gz"); writer->SetInput(m_SignalGen.m_MaskImage); writer->Update(); } catch(...) { MITK_INFO << "No mask image saved."; } setlocale(LC_ALL, currLocale.c_str()); } template< class ParameterType > ParameterType mitk::FiberfoxParameters::ReadVal(boost::property_tree::ptree::value_type const& v, std::string tag, ParameterType defaultValue, bool essential) { try { return v.second.get(tag); } catch (...) { if (essential) { mitkThrow() << "Parameter file corrupted. Essential tag is missing: '" << tag << "'"; } if (tag!="artifacts.noisetype") { MITK_INFO << "Tag '" << tag << "' not found. Using default value '" << defaultValue << "'."; m_MissingTags += "\n- "; m_MissingTags += tag; } return defaultValue; } } void mitk::FiberfoxParameters::UpdateSignalModels() { for (mitk::DiffusionSignalModel<>* m : m_FiberModelList) { m->SetGradientList(m_SignalGen.m_GradientDirections); m->SetBvalue(m_SignalGen.m_Bvalue); } for (mitk::DiffusionSignalModel<>* m : m_NonFiberModelList) { m->SetGradientList(m_SignalGen.m_GradientDirections); m->SetBvalue(m_SignalGen.m_Bvalue); } } void mitk::FiberfoxParameters::SetNumWeightedVolumes(int numGradients) { m_SignalGen.SetNumWeightedVolumes(numGradients); UpdateSignalModels(); } void mitk::FiberfoxParameters::SetGradienDirections(mitk::SignalGenerationParameters::GradientListType gradientList) { m_SignalGen.SetGradienDirections(gradientList); UpdateSignalModels(); } void mitk::FiberfoxParameters::SetGradienDirections(mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradientList) { m_SignalGen.SetGradienDirections(gradientList); UpdateSignalModels(); } void mitk::FiberfoxParameters::SetBvalue(double Bvalue) { m_SignalGen.m_Bvalue = Bvalue; UpdateSignalModels(); } void mitk::FiberfoxParameters::GenerateGradientHalfShell() { m_SignalGen.GenerateGradientHalfShell(); UpdateSignalModels(); } -void mitk::FiberfoxParameters::LoadParameters(std::string filename) +void mitk::FiberfoxParameters::LoadParameters(std::string filename, bool fix_seed) { - srand(time(0)); + if (fix_seed) + srand(0); + else + srand(time(0)); m_MissingTags = ""; if(filename.empty()) { return; } const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, nullptr ); if ( locale.compare(currLocale)!=0 ) { try { setlocale(LC_ALL, locale.c_str()); } catch(...) { MITK_INFO << "Could not set locale " << locale; } } boost::property_tree::ptree parameterTree; boost::property_tree::xml_parser::read_xml( filename, parameterTree ); m_FiberModelList.clear(); m_NonFiberModelList.clear(); if (m_NoiseModel) { m_NoiseModel = nullptr; } BOOST_FOREACH( boost::property_tree::ptree::value_type const& v1, parameterTree.get_child("fiberfox") ) { if( v1.first == "fibers" ) { m_Misc.m_CheckRealTimeFibersBox = ReadVal(v1,"realtime", m_Misc.m_CheckRealTimeFibersBox); m_Misc.m_CheckAdvancedFiberOptionsBox = ReadVal(v1,"showadvanced", m_Misc.m_CheckAdvancedFiberOptionsBox); m_Misc.m_CheckConstantRadiusBox = ReadVal(v1,"constantradius", m_Misc.m_CheckConstantRadiusBox); m_Misc.m_CheckIncludeFiducialsBox = ReadVal(v1,"includeFiducials", m_Misc.m_CheckIncludeFiducialsBox); switch (ReadVal(v1,"distribution", 0)) { case 0: m_FiberGen.m_Distribution = FiberGenerationParameters::DISTRIBUTE_UNIFORM; break; case 1: m_FiberGen.m_Distribution = FiberGenerationParameters::DISTRIBUTE_GAUSSIAN; break; default: m_FiberGen.m_Distribution = FiberGenerationParameters::DISTRIBUTE_UNIFORM; } m_FiberGen.m_Variance = ReadVal(v1,"variance", m_FiberGen.m_Variance); m_FiberGen.m_Density = ReadVal(v1,"density", m_FiberGen.m_Density); m_FiberGen.m_Sampling = ReadVal(v1,"spline.sampling", m_FiberGen.m_Sampling); m_FiberGen.m_Tension = ReadVal(v1,"spline.tension", m_FiberGen.m_Tension); m_FiberGen.m_Continuity = ReadVal(v1,"spline.continuity", m_FiberGen.m_Continuity); m_FiberGen.m_Bias = ReadVal(v1,"spline.bias", m_FiberGen.m_Bias); m_FiberGen.m_Rotation[0] = ReadVal(v1,"rotation.x", m_FiberGen.m_Rotation[0]); m_FiberGen.m_Rotation[1] = ReadVal(v1,"rotation.y", m_FiberGen.m_Rotation[1]); m_FiberGen.m_Rotation[2] = ReadVal(v1,"rotation.z", m_FiberGen.m_Rotation[2]); m_FiberGen.m_Translation[0] = ReadVal(v1,"translation.x", m_FiberGen.m_Translation[0]); m_FiberGen.m_Translation[1] = ReadVal(v1,"translation.y", m_FiberGen.m_Translation[1]); m_FiberGen.m_Translation[2] = ReadVal(v1,"translation.z", m_FiberGen.m_Translation[2]); m_FiberGen.m_Scale[0] = ReadVal(v1,"scale.x", m_FiberGen.m_Scale[0]); m_FiberGen.m_Scale[1] = ReadVal(v1,"scale.y", m_FiberGen.m_Scale[1]); m_FiberGen.m_Scale[2] = ReadVal(v1,"scale.z", m_FiberGen.m_Scale[2]); } else if ( v1.first == "image" ) { m_Misc.m_SignalModelString = ReadVal(v1,"signalmodelstring", m_Misc.m_SignalModelString); m_Misc.m_ArtifactModelString = ReadVal(v1,"artifactmodelstring", m_Misc.m_ArtifactModelString); m_Misc.m_OutputPath = ReadVal(v1,"outpath", m_Misc.m_OutputPath); m_Misc.m_CheckOutputVolumeFractionsBox = ReadVal(v1,"outputvolumefractions", m_Misc.m_CheckOutputVolumeFractionsBox); m_Misc.m_CheckAdvancedSignalOptionsBox = ReadVal(v1,"showadvanced", m_Misc.m_CheckAdvancedSignalOptionsBox); m_Misc.m_DoAddDistortions = ReadVal(v1,"artifacts.doAddDistortions", m_Misc.m_DoAddDistortions); m_Misc.m_DoAddNoise = ReadVal(v1,"artifacts.addnoise", m_Misc.m_DoAddNoise); m_Misc.m_DoAddGhosts = ReadVal(v1,"artifacts.addghosts", m_Misc.m_DoAddGhosts); m_Misc.m_DoAddAliasing = ReadVal(v1,"artifacts.addaliasing", m_Misc.m_DoAddAliasing); m_Misc.m_DoAddSpikes = ReadVal(v1,"artifacts.addspikes", m_Misc.m_DoAddSpikes); m_Misc.m_DoAddEddyCurrents = ReadVal(v1,"artifacts.addeddycurrents", m_Misc.m_DoAddEddyCurrents); m_SignalGen.m_ImageRegion.SetSize(0, ReadVal(v1,"basic.size.x",m_SignalGen.m_ImageRegion.GetSize(0))); m_SignalGen.m_ImageRegion.SetSize(1, ReadVal(v1,"basic.size.y",m_SignalGen.m_ImageRegion.GetSize(1))); m_SignalGen.m_ImageRegion.SetSize(2, ReadVal(v1,"basic.size.z",m_SignalGen.m_ImageRegion.GetSize(2))); m_SignalGen.m_ImageSpacing[0] = ReadVal(v1,"basic.spacing.x",m_SignalGen.m_ImageSpacing[0]); m_SignalGen.m_ImageSpacing[1] = ReadVal(v1,"basic.spacing.y",m_SignalGen.m_ImageSpacing[1]); m_SignalGen.m_ImageSpacing[2] = ReadVal(v1,"basic.spacing.z",m_SignalGen.m_ImageSpacing[2]); m_SignalGen.m_ImageOrigin[0] = ReadVal(v1,"basic.origin.x",m_SignalGen.m_ImageOrigin[0]); m_SignalGen.m_ImageOrigin[1] = ReadVal(v1,"basic.origin.y",m_SignalGen.m_ImageOrigin[1]); m_SignalGen.m_ImageOrigin[2] = ReadVal(v1,"basic.origin.z",m_SignalGen.m_ImageOrigin[2]); int i = 0; int j = 0; for(auto v : v1.second.get_child("basic.direction")) { m_SignalGen.m_ImageDirection[i][j] = boost::lexical_cast(v.second.data()); ++j; if (j==3) { j = 0; ++i; } } m_SignalGen.m_AcquisitionType = (SignalGenerationParameters::AcquisitionType)ReadVal(v1,"acquisitiontype", m_SignalGen.m_AcquisitionType); m_SignalGen.m_CoilSensitivityProfile = (SignalGenerationParameters::CoilSensitivityProfile)ReadVal(v1,"coilsensitivityprofile", m_SignalGen.m_CoilSensitivityProfile); m_SignalGen.m_NumberOfCoils = ReadVal(v1,"numberofcoils", m_SignalGen.m_NumberOfCoils); m_SignalGen.m_ReversePhase = ReadVal(v1,"reversephase", m_SignalGen.m_ReversePhase); m_SignalGen.m_PartialFourier = ReadVal(v1,"partialfourier", m_SignalGen.m_PartialFourier); m_SignalGen.m_NoiseVariance = ReadVal(v1,"noisevariance", m_SignalGen.m_NoiseVariance); m_SignalGen.m_tRep = ReadVal(v1,"trep", m_SignalGen.m_tRep); m_SignalGen.m_SignalScale = ReadVal(v1,"signalScale", m_SignalGen.m_SignalScale); m_SignalGen.m_tEcho = ReadVal(v1,"tEcho", m_SignalGen.m_tEcho); m_SignalGen.m_tLine = ReadVal(v1,"tLine", m_SignalGen.m_tLine); m_SignalGen.m_tInhom = ReadVal(v1,"tInhom", m_SignalGen.m_tInhom); m_SignalGen.m_SimulateKspaceAcquisition = ReadVal(v1,"simulatekspace", m_SignalGen.m_SimulateKspaceAcquisition); m_SignalGen.m_AxonRadius = ReadVal(v1,"axonRadius", m_SignalGen.m_AxonRadius); m_SignalGen.m_Spikes = ReadVal(v1,"artifacts.spikesnum", m_SignalGen.m_Spikes); m_SignalGen.m_SpikeAmplitude = ReadVal(v1,"artifacts.spikesscale", m_SignalGen.m_SpikeAmplitude); m_SignalGen.m_KspaceLineOffset = ReadVal(v1,"artifacts.kspaceLineOffset", m_SignalGen.m_KspaceLineOffset); m_SignalGen.m_EddyStrength = ReadVal(v1,"artifacts.eddyStrength", m_SignalGen.m_EddyStrength); m_SignalGen.m_Tau = ReadVal(v1,"artifacts.eddyTau", m_SignalGen.m_Tau); m_SignalGen.m_CroppingFactor = ReadVal(v1,"artifacts.aliasingfactor", m_SignalGen.m_CroppingFactor); m_SignalGen.m_Drift = ReadVal(v1,"artifacts.drift", m_SignalGen.m_Drift); m_SignalGen.m_DoAddGibbsRinging = ReadVal(v1,"artifacts.addringing", m_SignalGen.m_DoAddGibbsRinging); m_SignalGen.m_DoSimulateRelaxation = ReadVal(v1,"doSimulateRelaxation", m_SignalGen.m_DoSimulateRelaxation); m_SignalGen.m_DoDisablePartialVolume = ReadVal(v1,"doDisablePartialVolume", m_SignalGen.m_DoDisablePartialVolume); m_SignalGen.m_DoAddMotion = ReadVal(v1,"artifacts.doAddMotion", m_SignalGen.m_DoAddMotion); m_SignalGen.m_DoRandomizeMotion = ReadVal(v1,"artifacts.randomMotion", m_SignalGen.m_DoRandomizeMotion); m_SignalGen.m_DoAddDrift = ReadVal(v1,"artifacts.doAddDrift", m_SignalGen.m_DoAddDrift); m_SignalGen.m_Translation[0] = ReadVal(v1,"artifacts.translation0", m_SignalGen.m_Translation[0]); m_SignalGen.m_Translation[1] = ReadVal(v1,"artifacts.translation1", m_SignalGen.m_Translation[1]); m_SignalGen.m_Translation[2] = ReadVal(v1,"artifacts.translation2", m_SignalGen.m_Translation[2]); m_SignalGen.m_Rotation[0] = ReadVal(v1,"artifacts.rotation0", m_SignalGen.m_Rotation[0]); m_SignalGen.m_Rotation[1] = ReadVal(v1,"artifacts.rotation1", m_SignalGen.m_Rotation[1]); m_SignalGen.m_Rotation[2] = ReadVal(v1,"artifacts.rotation2", m_SignalGen.m_Rotation[2]); if (itksys::SystemTools::FileExists(filename+".bvals") && itksys::SystemTools::FileExists(filename+".bvecs")) { m_Misc.m_BvalsFile = filename+".bvals"; m_Misc.m_BvecsFile = filename+".bvecs"; m_SignalGen.SetGradienDirections( mitk::gradients::ReadBvalsBvecs(m_Misc.m_BvalsFile, m_Misc.m_BvecsFile, m_SignalGen.m_Bvalue) ); } else { m_SignalGen.m_Bvalue = ReadVal(v1,"bvalue", m_SignalGen.m_Bvalue); SignalGenerationParameters::GradientListType gradients; try { BOOST_FOREACH( boost::property_tree::ptree::value_type const& v2, v1.second.get_child("gradients") ) { SignalGenerationParameters::GradientType g; g[0] = ReadVal(v2,"x",0); g[1] = ReadVal(v2,"y",0); g[2] = ReadVal(v2,"z",0); gradients.push_back(g); } } catch(...) { MITK_INFO << "WARNING: Fiberfox parameters without any gradient directions loaded."; } m_SignalGen.SetGradienDirections(gradients); } m_Misc.m_MotionVolumesBox = ReadVal(v1,"artifacts.motionvolumes", m_Misc.m_MotionVolumesBox); m_SignalGen.m_MotionVolumes.clear(); if ( m_Misc.m_MotionVolumesBox == "random" ) { m_SignalGen.m_MotionVolumes.push_back(0); for ( size_t i=1; i < m_SignalGen.GetNumVolumes(); ++i ) { m_SignalGen.m_MotionVolumes.push_back( bool( rand()%2 ) ); } MITK_DEBUG << "mitkFiberfoxParameters.cpp: Case m_Misc.m_MotionVolumesBox == \"random\"."; } else if ( ! m_Misc.m_MotionVolumesBox.empty() ) { std::stringstream stream( m_Misc.m_MotionVolumesBox ); std::vector numbers; int nummer = std::numeric_limits::max(); while( stream >> nummer ) { if( nummer < std::numeric_limits::max() ) { numbers.push_back( nummer ); } } // If a list of negative numbers is given: if( *(std::min_element( numbers.begin(), numbers.end() )) < 0 && *(std::max_element( numbers.begin(), numbers.end() )) <= 0 ) // cave: -0 == +0 { for ( size_t i=0; i(m_SignalGen.GetNumVolumes()) && -number >= 0 ) m_SignalGen.m_MotionVolumes.at(-number) = false; } MITK_DEBUG << "mitkFiberfoxParameters.cpp: Case list of negative numbers."; } // If a list of positive numbers is given: else if( *(std::min_element( numbers.begin(), numbers.end() )) >= 0 && *(std::max_element( numbers.begin(), numbers.end() )) >= 0 ) { for ( size_t i=0; i(m_SignalGen.GetNumVolumes()) && number >= 0) m_SignalGen.m_MotionVolumes.at(number) = true; } MITK_DEBUG << "mitkFiberfoxParameters.cpp: Case list of positive numbers."; } else { MITK_WARN << "mitkFiberfoxParameters.cpp: Inconsistent list of numbers in m_MotionVolumesBox."; break; } } else { MITK_WARN << "mitkFiberfoxParameters.cpp: Cannot make sense of string in m_MotionVolumesBox."; break; } try { if (ReadVal(v1,"artifacts.noisetype","")=="rice") { m_NoiseModel = std::make_shared< mitk::RicianNoiseModel<> >(); m_NoiseModel->SetNoiseVariance(ReadVal(v1,"artifacts.noisevariance",m_NoiseModel->GetNoiseVariance())); } } catch(...) { MITK_DEBUG << "mitkFiberfoxParameters.cpp: caught some error while trying m_NoiseModel->SetNoiseVariance()"; // throw; } try { if (ReadVal(v1,"artifacts.noisetype","")=="chisquare") { m_NoiseModel = std::make_shared< mitk::ChiSquareNoiseModel<> >(); m_NoiseModel->SetNoiseVariance(ReadVal(v1,"artifacts.noisevariance",m_NoiseModel->GetNoiseVariance())); } } catch(...) { MITK_DEBUG << "mitkFiberfoxParameters.cpp: caught some error while trying m_NoiseModel->SetNoiseVariance()"; // throw; } BOOST_FOREACH( boost::property_tree::ptree::value_type const& v2, v1.second.get_child("compartments") ) { mitk::DiffusionSignalModel<>* signalModel = nullptr; std::string model = ReadVal(v2,"model","",true); if (model=="stick") { mitk::StickModel<>* model = new mitk::StickModel<>(); model->SetDiffusivity(ReadVal(v2,"d",model->GetDiffusivity())); model->SetT2(ReadVal(v2,"t2",model->GetT2())); model->SetT1(ReadVal(v2,"t1",model->GetT1())); model->SetBvalue(m_SignalGen.m_Bvalue); model->m_CompartmentId = ReadVal(v2,"ID",0,true); if (ReadVal(v2,"type","",true)=="fiber") m_FiberModelList.push_back(model); else if (ReadVal(v2,"type","",true)=="non-fiber") m_NonFiberModelList.push_back(model); signalModel = model; } else if (model=="tensor") { mitk::TensorModel<>* model = new mitk::TensorModel<>(); model->SetDiffusivity1(ReadVal(v2,"d1",model->GetDiffusivity1())); model->SetDiffusivity2(ReadVal(v2,"d2",model->GetDiffusivity2())); model->SetDiffusivity3(ReadVal(v2,"d3",model->GetDiffusivity3())); model->SetT2(ReadVal(v2,"t2",model->GetT2())); model->SetT1(ReadVal(v2,"t1",model->GetT1())); model->SetBvalue(m_SignalGen.m_Bvalue); model->m_CompartmentId = ReadVal(v2,"ID",0,true); if (ReadVal(v2,"type","",true)=="fiber") m_FiberModelList.push_back(model); else if (ReadVal(v2,"type","",true)=="non-fiber") m_NonFiberModelList.push_back(model); signalModel = model; } else if (model=="ball") { mitk::BallModel<>* model = new mitk::BallModel<>(); model->SetDiffusivity(ReadVal(v2,"d",model->GetDiffusivity())); model->SetT2(ReadVal(v2,"t2",model->GetT2())); model->SetT1(ReadVal(v2,"t1",model->GetT1())); model->SetBvalue(m_SignalGen.m_Bvalue); model->m_CompartmentId = ReadVal(v2,"ID",0,true); if (ReadVal(v2,"type","",true)=="fiber") m_FiberModelList.push_back(model); else if (ReadVal(v2,"type","",true)=="non-fiber") m_NonFiberModelList.push_back(model); signalModel = model; } else if (model=="astrosticks") { mitk::AstroStickModel<>* model = new AstroStickModel<>(); model->SetDiffusivity(ReadVal(v2,"d",model->GetDiffusivity())); model->SetT2(ReadVal(v2,"t2",model->GetT2())); model->SetT1(ReadVal(v2,"t1",model->GetT1())); model->SetBvalue(m_SignalGen.m_Bvalue); model->SetRandomizeSticks(ReadVal(v2,"randomize",model->GetRandomizeSticks())); model->m_CompartmentId = ReadVal(v2,"ID",0,true); if (ReadVal(v2,"type","",true)=="fiber") m_FiberModelList.push_back(model); else if (ReadVal(v2,"type","",true)=="non-fiber") m_NonFiberModelList.push_back(model); signalModel = model; } else if (model=="dot") { mitk::DotModel<>* model = new mitk::DotModel<>(); model->SetT2(ReadVal(v2,"t2",model->GetT2())); model->SetT1(ReadVal(v2,"t1",model->GetT1())); model->m_CompartmentId = ReadVal(v2,"ID",0,true); if (ReadVal(v2,"type","",true)=="fiber") m_FiberModelList.push_back(model); else if (ReadVal(v2,"type","",true)=="non-fiber") m_NonFiberModelList.push_back(model); signalModel = model; } else if (model=="prototype") { mitk::RawShModel<>* model = new mitk::RawShModel<>(); model->SetMaxNumKernels(ReadVal(v2,"maxNumSamples",model->GetMaxNumKernels())); model->SetFaRange(ReadVal(v2,"minFA",model->GetFaRange().first), ReadVal(v2,"maxFA",model->GetFaRange().second)); model->SetAdcRange(ReadVal(v2,"minADC",model->GetAdcRange().first), ReadVal(v2,"maxADC",model->GetAdcRange().second)); model->m_CompartmentId = ReadVal(v2,"ID",0,true); unsigned int numCoeffs = ReadVal(v2,"numCoeffs",0,true); unsigned int numSamples = ReadVal(v2,"numSamples",0,true); for (unsigned int j=0; j coeffs(numCoeffs); for (unsigned int k=0; k(v2,"kernels."+boost::lexical_cast(j)+".coeffs."+boost::lexical_cast(k),0,true); } model->SetShCoefficients( coeffs, ReadVal(v2,"kernels."+boost::lexical_cast(j)+".B0",0,true) ); } if (ReadVal(v2,"type","",true)=="fiber") { m_FiberModelList.push_back(model); } else if (ReadVal(v2,"type","",true)=="non-fiber") { m_NonFiberModelList.push_back(model); } // else ? signalModel = model; } if (signalModel!=nullptr) { try { itk::ImageFileReader::Pointer reader = itk::ImageFileReader::New(); if ( itksys::SystemTools::FileExists(filename+"_VOLUME"+ReadVal(v2,"ID","")+".nii.gz") ) reader->SetFileName(filename+"_VOLUME"+ReadVal(v2,"ID","")+".nii.gz"); else if ( itksys::SystemTools::FileExists(filename+"_VOLUME"+ReadVal(v2,"ID","")+".nii") ) reader->SetFileName(filename+"_VOLUME"+ReadVal(v2,"ID","")+".nii"); else reader->SetFileName(filename+"_VOLUME"+ReadVal(v2,"ID","")+".nrrd"); reader->Update(); signalModel->SetVolumeFractionImage(reader->GetOutput()); MITK_INFO << "Volume fraction image loaded for compartment " << signalModel->m_CompartmentId; } catch(...) { MITK_INFO << "No volume fraction image found for compartment " << signalModel->m_CompartmentId; } } } } else { } } UpdateSignalModels(); try { itk::ImageFileReader::Pointer reader = itk::ImageFileReader::New(); reader->SetFileName(filename+"_FMAP.nrrd"); if ( itksys::SystemTools::FileExists(filename+"_FMAP.nii.gz") ) reader->SetFileName(filename+"_FMAP.nii.gz"); else if ( itksys::SystemTools::FileExists(filename+"_FMAP.nii") ) reader->SetFileName(filename+"_FMAP.nii"); else reader->SetFileName(filename+"_FMAP.nrrd"); reader->Update(); m_SignalGen.m_FrequencyMap = reader->GetOutput(); MITK_INFO << "Frequency map loaded."; } catch(...) { MITK_INFO << "No frequency map found."; } try { itk::ImageFileReader::Pointer reader = itk::ImageFileReader::New(); if ( itksys::SystemTools::FileExists(filename+"_MASK.nii.gz") ) reader->SetFileName(filename+"_MASK.nii.gz"); else if ( itksys::SystemTools::FileExists(filename+"_MASK.nii") ) reader->SetFileName(filename+"_MASK.nii"); else reader->SetFileName(filename+"_MASK.nrrd"); reader->Update(); m_SignalGen.m_MaskImage = reader->GetOutput(); m_SignalGen.m_ImageRegion = m_SignalGen.m_MaskImage->GetLargestPossibleRegion(); m_SignalGen.m_ImageSpacing = m_SignalGen.m_MaskImage->GetSpacing(); m_SignalGen.m_ImageOrigin = m_SignalGen.m_MaskImage->GetOrigin(); m_SignalGen.m_ImageDirection = m_SignalGen.m_MaskImage->GetDirection(); MITK_INFO << "Mask image loaded."; } catch(...) { MITK_INFO << "No mask image found."; } setlocale(LC_ALL, currLocale.c_str()); } void mitk::FiberfoxParameters::PrintSelf() { MITK_INFO << "Not implemented :("; } diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h index 5ea3d0b7f0..5c495c6dba 100644 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h @@ -1,328 +1,328 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _MITK_FiberfoxParameters_H #define _MITK_FiberfoxParameters_H #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include namespace mitk { class MITKFIBERTRACKING_EXPORT FiberfoxParameters; /** Signal generation */ class MITKFIBERTRACKING_EXPORT SignalGenerationParameters { friend FiberfoxParameters; public: typedef itk::Image ItkFloatImgType; typedef itk::Image ItkUcharImgType; typedef itk::Vector GradientType; typedef std::vector GradientListType; enum CoilSensitivityProfile : int { COIL_CONSTANT, COIL_LINEAR, COIL_EXPONENTIAL }; enum AcquisitionType : int { SingleShotEpi, SpinEcho }; SignalGenerationParameters() : m_AcquisitionType(SignalGenerationParameters::SingleShotEpi) , m_SignalScale(100) , m_tEcho(100) , m_tRep(4000) , m_tLine(1) , m_tInhom(50) , m_ReversePhase(false) , m_PartialFourier(1.0) , m_NoiseVariance(0.001) , m_NumberOfCoils(1) , m_CoilSensitivityProfile(SignalGenerationParameters::COIL_CONSTANT) , m_SimulateKspaceAcquisition(false) , m_AxonRadius(0) , m_DoDisablePartialVolume(false) , m_Spikes(0) , m_SpikeAmplitude(1) , m_KspaceLineOffset(0) , m_EddyStrength(300) , m_Tau(70) , m_CroppingFactor(1) , m_Drift(0.06) , m_DoAddGibbsRinging(false) , m_DoSimulateRelaxation(true) , m_DoAddMotion(false) , m_DoRandomizeMotion(true) , m_DoAddDrift(false) , m_FrequencyMap(nullptr) , m_MaskImage(nullptr) , m_Bvalue(1000) { m_ImageRegion.SetSize(0, 12); m_ImageRegion.SetSize(1, 12); m_ImageRegion.SetSize(2, 3); m_ImageSpacing.Fill(2.0); m_ImageOrigin.Fill(0.0); m_ImageDirection.SetIdentity(); m_Translation.Fill(0.0); m_Rotation.Fill(0.0); SetNumWeightedVolumes(6); } /** input/output image specifications */ itk::ImageRegion<3> m_CroppedRegion; ///< Image size with reduced FOV. itk::ImageRegion<3> m_ImageRegion; ///< Image size. itk::Vector m_ImageSpacing; ///< Image voxel size. itk::Point m_ImageOrigin; ///< Image origin. itk::Matrix m_ImageDirection; ///< Image rotation matrix. /** Other acquisitions parameters */ AcquisitionType m_AcquisitionType; ///< determines k-space trajectory and maximum echo position(s) float m_SignalScale; ///< Scaling factor for output signal (before noise is added). float m_tEcho; ///< Echo time TE. float m_tRep; ///< Echo time TR. float m_tLine; ///< k-space line readout time (dwell time). float m_tInhom; ///< T2' bool m_ReversePhase; ///< If true, the phase readout direction will be inverted (-y instead of y) float m_PartialFourier; ///< Partial fourier factor (0.5-1) float m_NoiseVariance; ///< Variance of complex gaussian noise int m_NumberOfCoils; ///< Number of coils in multi-coil acquisition CoilSensitivityProfile m_CoilSensitivityProfile; ///< Choose between constant, linear or exponential sensitivity profile of the used coils bool m_SimulateKspaceAcquisition;///< Flag to enable/disable k-space acquisition simulation double m_AxonRadius; ///< Determines compartment volume fractions (0 == automatic axon radius estimation) bool m_DoDisablePartialVolume; ///< Disable partial volume effects. Each voxel is either all fiber or all non-fiber. /** Artifacts and other effects */ unsigned int m_Spikes; ///< Number of spikes randomly appearing in the image float m_SpikeAmplitude; ///< amplitude of spikes relative to the largest signal intensity (magnitude of complex) float m_KspaceLineOffset; ///< Causes N/2 ghosts. Larger offset means stronger ghost. float m_EddyStrength; ///< Strength of eddy current induced gradients in mT/m. float m_Tau; ///< Eddy current decay constant (in ms) float m_CroppingFactor; ///< FOV size in y-direction is multiplied by this factor. Causes aliasing artifacts. float m_Drift; ///< Global signal decrease by the end of the acquisition. bool m_DoAddGibbsRinging; ///< Add Gibbs ringing artifact bool m_DoSimulateRelaxation; ///< Add T2 relaxation effects bool m_DoAddMotion; ///< Enable motion artifacts. bool m_DoRandomizeMotion; ///< Toggles between random and linear motion. bool m_DoAddDrift; ///< Add quadratic signal drift. std::vector< bool > m_MotionVolumes; ///< Indicates the image volumes that are affected by motion ///< with positive numbers, inverted logic with negative numbers. itk::Vector m_Translation; ///< Maximum translational motion. itk::Vector m_Rotation; ///< Maximum rotational motion. ItkFloatImgType::Pointer m_FrequencyMap; ///< If != nullptr, distortions are added to the image using this frequency map. ItkUcharImgType::Pointer m_MaskImage; ///< Signal is only genrated inside of the mask image. std::vector< int > GetBaselineIndices(); ///< Returns list of nun-diffusion-weighted image volume indices unsigned int GetFirstBaselineIndex(); ///< Returns index of first non-diffusion-weighted image volume bool IsBaselineIndex(unsigned int idx); ///< Checks if image volume with given index is non-diffusion-weighted volume or not. unsigned int GetNumWeightedVolumes(); ///< Get number of diffusion-weighted image volumes unsigned int GetNumBaselineVolumes(); ///< Get number of non-diffusion-weighted image volumes unsigned int GetNumVolumes(); ///< Get number of baseline and diffusion-weighted image volumes GradientListType GetGradientDirections(); ///< Return gradient direction container mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer GetItkGradientContainer(); GradientType GetGradientDirection(unsigned int i); std::vector< int > GetBvalues(); ///< Returns a vector with all unique b-values (determined by the gradient magnitudes) double GetBvalue(); void ApplyDirectionMatrix(); protected: unsigned int m_NumGradients; ///< Number of diffusion-weighted image volumes. unsigned int m_NumBaseline; ///< Number of non-diffusion-weighted image volumes. GradientListType m_GradientDirections; ///< Total number of image volumes. double m_Bvalue; ///< Acquisition b-value void SetNumWeightedVolumes(int numGradients); ///< Automaticall calls GenerateGradientHalfShell() afterwards. void SetGradienDirections(GradientListType gradientList); void SetGradienDirections(mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradientList); void GenerateGradientHalfShell(); ///< Generates half shell of gradient directions (with m_NumGradients non-zero directions) }; /** Fiber generation */ class MITKFIBERTRACKING_EXPORT FiberGenerationParameters { public: enum FiberDistribution { DISTRIBUTE_UNIFORM, // distribute fibers uniformly in the ROIs DISTRIBUTE_GAUSSIAN // distribute fibers using a 2D gaussian }; typedef std::vector< std::vector< mitk::PlanarEllipse::Pointer > > FiducialListType; typedef std::vector< std::vector< unsigned int > > FlipListType; FiberGenerationParameters() : m_Distribution(DISTRIBUTE_UNIFORM) , m_Density(100) , m_Variance(100) , m_Sampling(1) , m_Tension(0) , m_Continuity(0) , m_Bias(0) { m_Rotation.Fill(0.0); m_Translation.Fill(0.0); m_Scale.Fill(1.0); } FiberDistribution m_Distribution; unsigned int m_Density; double m_Variance; double m_Sampling; double m_Tension; double m_Continuity; double m_Bias; mitk::Vector3D m_Rotation; mitk::Vector3D m_Translation; mitk::Vector3D m_Scale; FlipListType m_FlipList; ///< contains flags indicating a flip of the 2D fiber x-coordinates (needed to resolve some unwanted fiber twisting) FiducialListType m_Fiducials; ///< container of the planar ellipses used as fiducials for the fiber generation process }; /** GUI persistence, input, output, ... */ class MITKFIBERTRACKING_EXPORT MiscFiberfoxParameters { public: MiscFiberfoxParameters() : m_ResultNode(DataNode::New()) , m_ParentNode(nullptr) , m_SignalModelString("") , m_ArtifactModelString("") , m_OutputPath("/tmp/") , m_OutputPrefix("fiberfox") , m_AfterSimulationMessage("") , m_BvalsFile("") , m_BvecsFile("") , m_CheckOutputVolumeFractionsBox(false) , m_CheckAdvancedSignalOptionsBox(false) , m_DoAddNoise(false) , m_DoAddGhosts(false) , m_DoAddAliasing(false) , m_DoAddSpikes(false) , m_DoAddEddyCurrents(false) , m_DoAddDistortions(false) , m_MotionVolumesBox("random") , m_CheckRealTimeFibersBox(true) , m_CheckAdvancedFiberOptionsBox(false) , m_CheckConstantRadiusBox(false) , m_CheckIncludeFiducialsBox(true) {} DataNode::Pointer m_ResultNode; ///< Stores resulting image. DataNode::Pointer m_ParentNode; ///< Parent node of result node. std::string m_SignalModelString; ///< Appendet to the name of the result node std::string m_ArtifactModelString; ///< Appendet to the name of the result node std::string m_OutputPath; ///< Image is automatically saved to the specified folder after simulation is finished. std::string m_OutputPrefix; /** Prefix for filename of output files and logfile. */ std::string m_AfterSimulationMessage; ///< Store messages that are displayed after the simulation has finished (e.g. warnings, automatic parameter adjustments etc.) std::string m_BvalsFile; std::string m_BvecsFile; /** member variables that store the check-state of GUI checkboxes */ // image generation bool m_CheckOutputVolumeFractionsBox; bool m_CheckAdvancedSignalOptionsBox; bool m_DoAddNoise; bool m_DoAddGhosts; bool m_DoAddAliasing; bool m_DoAddSpikes; bool m_DoAddEddyCurrents; bool m_DoAddDistortions; std::string m_MotionVolumesBox; // fiber generation bool m_CheckRealTimeFibersBox; bool m_CheckAdvancedFiberOptionsBox; bool m_CheckConstantRadiusBox; bool m_CheckIncludeFiducialsBox; }; /** * \brief Datastructure to manage the Fiberfox signal generation parameters. * */ class MITKFIBERTRACKING_EXPORT FiberfoxParameters { public: typedef itk::Image ItkFloatImgType; typedef itk::Image ItkDoubleImgType; typedef itk::Image ItkUcharImgType; typedef DiffusionSignalModel DiffusionModelType; typedef std::vector< DiffusionModelType* > DiffusionModelListType; typedef DiffusionNoiseModel NoiseModelType; FiberfoxParameters(); FiberfoxParameters(const FiberfoxParameters ¶ms); ~FiberfoxParameters(); /** Not templated parameters */ FiberGenerationParameters m_FiberGen; ///< Fiber generation parameters SignalGenerationParameters m_SignalGen; ///< Signal generation parameters MiscFiberfoxParameters m_Misc; ///< GUI realted and I/O parameters /** Templated parameters */ DiffusionModelListType m_FiberModelList; ///< Intra- and inter-axonal compartments. DiffusionModelListType m_NonFiberModelList; ///< Extra-axonal compartments. std::shared_ptr< NoiseModelType > m_NoiseModel; ///< If != nullptr, noise is added to the image. void GenerateGradientHalfShell(); void SetNumWeightedVolumes(int numGradients); ///< Automaticall calls GenerateGradientHalfShell() afterwards. void SetGradienDirections(mitk::SignalGenerationParameters::GradientListType gradientList); void SetGradienDirections(mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradientList); void SetBvalue(double Bvalue); void UpdateSignalModels(); void ClearFiberParameters(); void ClearSignalParameters(); void ApplyDirectionMatrix(); void PrintSelf(); ///< Print parameters to stdout. void SaveParameters(std::string filename); ///< Save image generation parameters to .ffp file. - void LoadParameters(std::string filename); ///< Load image generation parameters from .ffp file. + void LoadParameters(std::string filename, bool fix_seed=false); ///< Load image generation parameters from .ffp file. template< class ParameterType > ParameterType ReadVal(boost::property_tree::ptree::value_type const& v, std::string tag, ParameterType defaultValue, bool essential=false); std::string m_MissingTags; }; } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Testing/CMakeLists.txt b/Modules/DiffusionImaging/FiberTracking/Testing/CMakeLists.txt index b308270069..28d176e5d9 100644 --- a/Modules/DiffusionImaging/FiberTracking/Testing/CMakeLists.txt +++ b/Modules/DiffusionImaging/FiberTracking/Testing/CMakeLists.txt @@ -1,24 +1,25 @@ MITK_CREATE_MODULE_TESTS() if("${CMAKE_SIZEOF_VOID_P}" EQUAL "8") mitkAddCustomModuleTest(mitkFiberBundleReaderWriterTest mitkFiberBundleReaderWriterTest) # Temporarily disabled. Since method relies on random numbers, the behaviour is not consistent across different systems. Solution? #mitkAddCustomModuleTest(mitkGibbsTrackingTest mitkGibbsTrackingTest ${MITK_DATA_DIR}/DiffusionImaging/qBallImage.qbi ${MITK_DATA_DIR}/DiffusionImaging/diffusionImageMask.nrrd ${MITK_DATA_DIR}/DiffusionImaging/gibbsTrackingParameters.gtp ${MITK_DATA_DIR}/DiffusionImaging/gibbsTractogram.fib) mitkAddCustomModuleTest(mitkFiberTransformationTest mitkFiberTransformationTest ${MITK_DATA_DIR}/DiffusionImaging/fiberBundleX.fib ${MITK_DATA_DIR}/DiffusionImaging/fiberBundleX_transformed.fib) mitkAddCustomModuleTest(mitkFiberExtractionTest mitkFiberExtractionTest ${MITK_DATA_DIR}/DiffusionImaging/fiberBundleX.fib ${MITK_DATA_DIR}/DiffusionImaging/fiberBundleX_extracted.fib ${MITK_DATA_DIR}/DiffusionImaging/ROI1.pf ${MITK_DATA_DIR}/DiffusionImaging/ROI2.pf ${MITK_DATA_DIR}/DiffusionImaging/ROI3.pf ${MITK_DATA_DIR}/DiffusionImaging/ROIIMAGE.nrrd ${MITK_DATA_DIR}/DiffusionImaging/fiberBundleX_inside.fib ${MITK_DATA_DIR}/DiffusionImaging/fiberBundleX_outside.fib ${MITK_DATA_DIR}/DiffusionImaging/fiberBundleX_passing-mask.fib ${MITK_DATA_DIR}/DiffusionImaging/fiberBundleX_ending-in-mask.fib ${MITK_DATA_DIR}/DiffusionImaging/fiberBundleX_subtracted.fib ${MITK_DATA_DIR}/DiffusionImaging/fiberBundleX_added.fib) mitkAddCustomModuleTest(mitkFiberGenerationTest mitkFiberGenerationTest ${MITK_DATA_DIR}/DiffusionImaging/Fiberfox/Fiducial_0.pf ${MITK_DATA_DIR}/DiffusionImaging/Fiberfox/Fiducial_1.pf ${MITK_DATA_DIR}/DiffusionImaging/Fiberfox/Fiducial_2.pf ${MITK_DATA_DIR}/DiffusionImaging/Fiberfox/uniform.fib ${MITK_DATA_DIR}/DiffusionImaging/Fiberfox/gaussian.fib) +mitkAddCustomModuleTest(mitkFiberfoxSignalGenerationBrainSliceTest mitkFiberfoxSignalGenerationBrainSliceTest) mitkAddCustomModuleTest(mitkFiberfoxSignalGenerationTest mitkFiberfoxSignalGenerationTest) mitkAddCustomModuleTest(mitkMachineLearningTrackingTest mitkMachineLearningTrackingTest) mitkAddCustomModuleTest(mitkStreamlineTractographyTest mitkStreamlineTractographyTest) mitkAddCustomModuleTest(mitkFiberProcessingTest mitkFiberProcessingTest) mitkAddCustomModuleTest(mitkFiberFitTest mitkFiberFitTest) mitkAddCustomModuleTest(mitkPeakShImageReaderTest mitkPeakShImageReaderTest) if(MITK_ENABLE_RENDERING_TESTING) # apparently does not work on ubuntu mitkAddCustomModuleTest(mitkFiberMapper3DTest mitkFiberMapper3DTest) ENDIF() ENDIF() diff --git a/Modules/DiffusionImaging/FiberTracking/Testing/files.cmake b/Modules/DiffusionImaging/FiberTracking/Testing/files.cmake index ff08d834ad..3d06aa5391 100644 --- a/Modules/DiffusionImaging/FiberTracking/Testing/files.cmake +++ b/Modules/DiffusionImaging/FiberTracking/Testing/files.cmake @@ -1,16 +1,17 @@ SET(MODULE_CUSTOM_TESTS mitkFiberBundleReaderWriterTest.cpp mitkGibbsTrackingTest.cpp mitkStreamlineTractographyTest.cpp mitkFiberTransformationTest.cpp mitkFiberExtractionTest.cpp mitkFiberGenerationTest.cpp mitkFiberfoxSignalGenerationTest.cpp + mitkFiberfoxSignalGenerationBrainSliceTest.cpp mitkMachineLearningTrackingTest.cpp mitkFiberProcessingTest.cpp mitkFiberFitTest.cpp mitkFiberMapper3DTest.cpp mitkPeakShImageReaderTest.cpp ) diff --git a/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxSignalGenerationBrainSliceTest.cpp b/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxSignalGenerationBrainSliceTest.cpp new file mode 100644 index 0000000000..4cbc7cc057 --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxSignalGenerationBrainSliceTest.cpp @@ -0,0 +1,205 @@ +/*=================================================================== + +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 +#include +#include +#include +#include +#include +#include + +#include "mitkTestFixture.h" + +class mitkFiberfoxSignalGenerationBrainSliceTestSuite : public mitk::TestFixture +{ + + CPPUNIT_TEST_SUITE(mitkFiberfoxSignalGenerationBrainSliceTestSuite); + MITK_TEST(Test0); + MITK_TEST(Test1); + MITK_TEST(Test2); + CPPUNIT_TEST_SUITE_END(); + + typedef itk::VectorImage< short, 3> ItkDwiType; + +private: + +public: + + /** Members used inside the different (sub-)tests. All members are initialized via setUp().*/ + FiberBundle::Pointer m_FiberBundle; + mitk::Image::Pointer m_Template; + std::vector< FiberfoxParameters > m_Parameters; + std::vector< mitk::Image::Pointer > m_RefImages; + + void setUp() override + { +// std::srand(0); +// omp_set_num_threads(1); +// mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); + + m_FiberBundle = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/SignalGen_BrainSlice.fib")); + m_Template = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/SignalGen_BrainSliceTemplate.nii.gz")); + + { + FiberfoxParameters parameters; + parameters.LoadParameters(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/BrainSlice_1.ffp"), true); + adjust_to_template(parameters); + m_Parameters.push_back(parameters); + m_RefImages.push_back(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/BrainSlice_1_OUT.nii.gz"))); + } + + { + FiberfoxParameters parameters; + parameters.LoadParameters(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/BrainSlice_2.ffp"), true); + adjust_to_template(parameters); + m_Parameters.push_back(parameters); + m_RefImages.push_back(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/BrainSlice_2_OUT.nii.gz"))); + } + + { + FiberfoxParameters parameters; + parameters.LoadParameters(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/BrainSlice_3.ffp"), true); + adjust_to_template(parameters); + m_Parameters.push_back(parameters); + m_RefImages.push_back(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/BrainSlice_3_OUT.nii.gz"))); + } + } + + void adjust_to_template(mitk::FiberfoxParameters& parameters) + { + typedef itk::Image< short, 3 > ItkImageType; + ItkImageType::Pointer itkImagePointer = ItkImageType::New(); + mitk::CastToItkImage(m_Template, itkImagePointer); + parameters.m_SignalGen.m_ImageRegion = itkImagePointer->GetLargestPossibleRegion(); + parameters.m_SignalGen.m_ImageSpacing = itkImagePointer->GetSpacing(); + parameters.m_SignalGen.m_ImageOrigin = itkImagePointer->GetOrigin(); + parameters.m_SignalGen.m_ImageDirection = itkImagePointer->GetDirection(); + } + + void tearDown() override + { + + } + + bool CompareDwi(itk::VectorImage< short, 3 >* dwi1, itk::VectorImage< short, 3 >* dwi2) + { + bool out = true; + typedef itk::VectorImage< short, 3 > DwiImageType; + try{ + itk::ImageRegionIterator< DwiImageType > it1(dwi1, dwi1->GetLargestPossibleRegion()); + itk::ImageRegionIterator< DwiImageType > it2(dwi2, dwi2->GetLargestPossibleRegion()); + int count = 0; + while(!it1.IsAtEnd()) + { + for (unsigned int i=0; iGetVectorLength(); ++i) + { + short d = abs(it1.Get()[i]-it2.Get()[i]); + + if (d>0) + { + if (count<10) + { + MITK_INFO << "**************************************"; + MITK_INFO << "Test value: " << it1.GetIndex() << ":" << it1.Get()[i]; + MITK_INFO << "Ref. value: " << it2.GetIndex() << ":" << it2.Get()[i]; + } + out = false; + count++; + } + } + ++it1; + ++it2; + } + if (count>=10) + MITK_INFO << "Skipping errors."; + MITK_INFO << "Errors detected: " << count; + } + catch(...) + { + return false; + } + return out; + } + + void StartSimulation(FiberfoxParameters parameters, mitk::Image::Pointer refImage, std::string out) + { + itk::TractsToDWIImageFilter< short >::Pointer tractsToDwiFilter = itk::TractsToDWIImageFilter< short >::New(); + tractsToDwiFilter->SetUseConstantRandSeed(true); + tractsToDwiFilter->SetParameters(parameters); + tractsToDwiFilter->SetFiberBundle(m_FiberBundle); + tractsToDwiFilter->Update(); + + mitk::Image::Pointer testImage = mitk::GrabItkImageMemory( tractsToDwiFilter->GetOutput() ); + mitk::DiffusionPropertyHelper::SetGradientContainer(testImage, parameters.m_SignalGen.GetItkGradientContainer()); + mitk::DiffusionPropertyHelper::SetReferenceBValue(testImage, parameters.m_SignalGen.GetBvalue()); + + mitk::DiffusionPropertyHelper::InitializeImage( testImage ); + + if (refImage.IsNotNull()) + { + if(mitk::DiffusionPropertyHelper::GetGradientContainer(refImage).IsNotNull() ) + { + ItkDwiType::Pointer itkTestImagePointer = ItkDwiType::New(); + mitk::CastToItkImage(testImage, itkTestImagePointer); + ItkDwiType::Pointer itkRefImagePointer = ItkDwiType::New(); + mitk::CastToItkImage(refImage, itkRefImagePointer); + + bool cond = CompareDwi(itkTestImagePointer, itkRefImagePointer); + if (!cond) + { + MITK_INFO << "Saving test image to " << mitk::IOUtil::GetTempPath(); + mitk::IOUtil::Save(testImage, mitk::IOUtil::GetTempPath()+out); + } + CPPUNIT_ASSERT_MESSAGE("Simulated images should be equal", cond); + } + } + } + + void Test0() + { + StartSimulation(m_Parameters.at(0), m_RefImages.at(0), "BrainSlice_1_OUT.dwi"); + } + + void Test1() + { + StartSimulation(m_Parameters.at(1), m_RefImages.at(1), "BrainSlice_2_OUT.dwi"); + } + + void Test2() + { + StartSimulation(m_Parameters.at(2), m_RefImages.at(2), "BrainSlice_3_OUT.dwi"); + } + +}; + +MITK_TEST_SUITE_REGISTRATION(mitkFiberfoxSignalGenerationBrainSlice)