diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxNoTemplateParameters.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxNoTemplateParameters.cpp new file mode 100644 index 0000000000..e3c89ef9ed --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxNoTemplateParameters.cpp @@ -0,0 +1,150 @@ +/*=================================================================== + +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 "mitkFiberfoxNoTemplateParameters.h" + +void mitk::SignalGenerationParameters::GenerateGradientHalfShell() +{ + int NPoints = 2*m_NumGradients; + m_GradientDirections.clear(); + + m_NumBaseline = NPoints/20; + if (m_NumBaseline==0) + m_NumBaseline=1; + + GradientType g; + g.Fill(0.0); + for (unsigned int i=0; i theta; theta.set_size(NPoints); + vnl_vector phi; phi.set_size(NPoints); + double C = sqrt(4*M_PI); + phi(0) = 0.0; + phi(NPoints-1) = 0.0; + for(int i=0; i0 && i 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::SignalGenerationParameters::GradientType mitk::SignalGenerationParameters::GetGradientDirection(unsigned int i) +{ + return m_GradientDirections.at(i); +} + +void mitk::SignalGenerationParameters::SetNumWeightedVolumes(int numGradients) +{ + m_NumGradients = numGradients; + GenerateGradientHalfShell(); +} + +void mitk::SignalGenerationParameters::SetGradienDirections(GradientListType gradientList) +{ + m_GradientDirections = gradientList; + m_NumGradients = 0; + m_NumBaseline = 0; + for( unsigned int i=0; im_GradientDirections.size(); i++) + { + if (m_GradientDirections.at(i).GetNorm()>0.0001) + m_NumGradients++; + else + m_NumBaseline++; + } +} + +void mitk::SignalGenerationParameters::SetGradienDirections(mitk::DiffusionImage::GradientDirectionContainerType::Pointer gradientList) +{ + m_NumGradients = 0; + m_NumBaseline = 0; + m_GradientDirections.clear(); + + for( unsigned int i=0; iSize(); i++) + { + GradientType g; + g[0] = gradientList->at(i)[0]; + g[1] = gradientList->at(i)[1]; + g[2] = gradientList->at(i)[2]; + m_GradientDirections.push_back(g); + + if (m_GradientDirections.at(i).GetNorm()>0.0001) + m_NumGradients++; + else + m_NumBaseline++; + } +} diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxNoTemplateParameters.h similarity index 57% copy from Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h copy to Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxNoTemplateParameters.h index 4e43573b4b..085e571b7f 100644 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxNoTemplateParameters.h @@ -1,321 +1,228 @@ /*=================================================================== 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 +#ifndef _MITK_FiberfoxNoTemplateParameters_H +#define _MITK_FiberfoxNoTemplateParameters_H #include -#include -#include -#include +#include #include -#include -#include #include #include -#include -#include -#include -#include -#include -#include - -using namespace std; - namespace mitk { /** Signal generation */ -class MitkFiberTracking_EXPORT SignalGenerationParameters +class SignalGenerationParameters { public: typedef itk::Image ItkDoubleImgType; typedef itk::Image ItkUcharImgType; typedef itk::Vector GradientType; typedef std::vector GradientListType; enum DiffusionDirectionMode { FIBER_TANGENT_DIRECTIONS, MAIN_FIBER_DIRECTIONS, RANDOM_DIRECTIONS }; SignalGenerationParameters() - : m_FrequencyMap(NULL) - , m_MaskImage(NULL) - , m_SignalScale(100) + : m_SignalScale(100) , m_tEcho(100) , m_tLine(1) , m_tInhom(50) , m_Bvalue(1000) , m_SimulateKspaceAcquisition(false) , m_AxonRadius(0) , m_DiffusionDirectionMode(SignalGenerationParameters::FIBER_TANGENT_DIRECTIONS) , m_FiberSeparationThreshold(30) , m_Spikes(0) , m_SpikeAmplitude(1) , m_KspaceLineOffset(0) , m_EddyStrength(0) , m_Tau(70) , m_CroppingFactor(1) , m_DoAddGibbsRinging(false) , m_DoSimulateRelaxation(true) , m_DoDisablePartialVolume(false) , m_DoAddMotion(false) , m_DoRandomizeMotion(true) + , m_FrequencyMap(NULL) + , m_MaskImage(NULL) { 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); } + + ~SignalGenerationParameters() + { + } + /** input/output image specifications */ itk::ImageRegion<3> m_ImageRegion; ///< Image size. itk::Vector m_ImageSpacing; ///< Image voxel size. itk::Point m_ImageOrigin; ///< Image origin. itk::Matrix m_ImageDirection; ///< Image rotation matrix. /** Other acquisitions parameters */ double m_SignalScale; ///< Scaling factor for output signal (before noise is added). double m_tEcho; ///< Echo time TE. double m_tLine; ///< k-space line readout time. double m_tInhom; ///< T2' double m_Bvalue; ///< Acquisition b-value bool m_SimulateKspaceAcquisition;///< 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. DiffusionDirectionMode m_DiffusionDirectionMode; ///< Determines how the main diffusion direction of the signal models is selected double m_FiberSeparationThreshold; ///< Used for random and and mein fiber deriction DiffusionDirectionMode /** Artifacts and other effects */ unsigned int m_Spikes; ///< Number of spikes randomly appearing in the image double m_SpikeAmplitude; ///< amplitude of spikes relative to the largest signal intensity (magnitude of complex) double m_KspaceLineOffset; ///< Causes N/2 ghosts. Larger offset means stronger ghost. double m_EddyStrength; ///< Strength of eddy current induced gradients in mT/m. double m_Tau; ///< Eddy current decay constant (in ms) double m_CroppingFactor; ///< FOV size in y-direction is multiplied by this factor. Causes aliasing artifacts. 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. itk::Vector m_Translation; ///< Maximum translational motion. itk::Vector m_Rotation; ///< Maximum rotational motion. ItkDoubleImgType::Pointer m_FrequencyMap; ///< If != NULL, distortions are added to the image using this frequency map. ItkUcharImgType::Pointer m_MaskImage; ///< Signal is only genrated inside of the mask image. - inline void GenerateGradientHalfShell(); ///< Generates half shell of gradient directions (with m_NumGradients non-zero directions) - inline std::vector< int > GetBaselineIndices(); ///< Returns list of nun-diffusion-weighted image volume indices - inline unsigned int GetFirstBaselineIndex(); ///< Returns index of first non-diffusion-weighted image volume - inline bool IsBaselineIndex(unsigned int idx); ///< Checks if image volume with given index is non-diffusion-weighted volume or not. - inline unsigned int GetNumWeightedVolumes(); ///< Get number of diffusion-weighted image volumes - inline unsigned int GetNumBaselineVolumes(); ///< Get number of non-diffusion-weighted image volumes - inline unsigned int GetNumVolumes(); ///< Get number of baseline and diffusion-weighted image volumes - inline GradientListType GetGradientDirections(); ///< Return gradient direction container - inline GradientType GetGradientDirection(unsigned int i); + void GenerateGradientHalfShell(); ///< Generates half shell of gradient directions (with m_NumGradients non-zero directions) + 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 + GradientType GetGradientDirection(unsigned int i); - inline void SetNumWeightedVolumes(int numGradients); ///< Automaticall calls GenerateGradientHalfShell() afterwards. - inline void SetGradienDirections(GradientListType gradientList); - inline void SetGradienDirections(mitk::DiffusionImage::GradientDirectionContainerType::Pointer gradientList); + void SetNumWeightedVolumes(int numGradients); ///< Automaticall calls GenerateGradientHalfShell() afterwards. + void SetGradienDirections(GradientListType gradientList); + void SetGradienDirections(mitk::DiffusionImage::GradientDirectionContainerType::Pointer gradientList); protected: unsigned int m_NumGradients; ///< Number of diffusion-weighted image volumes. unsigned int m_NumBaseline; ///< Number of non-diffusion-weighted image volumes. GradientListType m_GradientDirections; ///< Total number of image volumes. }; /** Fiber generation */ -class MitkFiberTracking_EXPORT FiberGenerationParameters +class FiberGenerationParameters { public: enum FiberDistribution{ DISTRIBUTE_UNIFORM, // distribute fibers uniformly in the ROIs DISTRIBUTE_GAUSSIAN // distribute fibers using a 2D gaussian }; - typedef vector< vector< mitk::PlanarEllipse::Pointer > > FiducialListType; - typedef vector< vector< unsigned int > > FlipListType; + 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 +class MiscFiberfoxParameters { public: MiscFiberfoxParameters() : m_ResultNode(DataNode::New()) , m_ParentNode(NULL) , m_SignalModelString("") , m_ArtifactModelString("") , m_OutputPath("") , m_CheckOutputVolumeFractionsBox(false) , m_CheckAdvancedSignalOptionsBox(false) , m_CheckAddNoiseBox(false) , m_CheckAddGhostsBox(false) , m_CheckAddAliasingBox(false) , m_CheckAddSpikesBox(false) , m_CheckAddEddyCurrentsBox(false) , m_CheckAddDistortionsBox(false) , 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. - string m_SignalModelString; ///< Appendet to the name of the result node - string m_ArtifactModelString; ///< Appendet to the name of the result node - string m_OutputPath; ///< Image is automatically saved to the specified folder after simulation is finished. + 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. /** member variables that store the check-state of GUI checkboxes */ // image generation bool m_CheckOutputVolumeFractionsBox; bool m_CheckAdvancedSignalOptionsBox; bool m_CheckAddNoiseBox; bool m_CheckAddGhostsBox; bool m_CheckAddAliasingBox; bool m_CheckAddSpikesBox; bool m_CheckAddEddyCurrentsBox; bool m_CheckAddDistortionsBox; // fiber generation bool m_CheckRealTimeFibersBox; bool m_CheckAdvancedFiberOptionsBox; bool m_CheckConstantRadiusBox; bool m_CheckIncludeFiducialsBox; }; - -/** - * \brief Datastructure to manage the Fiberfox signal generation parameters. - * - */ -template< class ScalarType = double > -class MitkFiberTracking_EXPORT FiberfoxParameters -{ -public: - - typedef itk::Image ItkDoubleImgType; - typedef itk::Image ItkUcharImgType; - typedef DiffusionSignalModel DiffusionModelType; - typedef std::vector< DiffusionModelType* > DiffusionModelListType; - typedef DiffusionNoiseModel NoiseModelType; - - FiberfoxParameters(); - ~FiberfoxParameters(); - - /** Get same parameter object with different template parameter */ - template< class OutType > - FiberfoxParameters< OutType > CopyParameters() - { - FiberfoxParameters< OutType > out; - - out.m_FiberGen = m_FiberGen; - out.m_SignalGen = m_SignalGen; - out.m_Misc = m_Misc; - - if (m_NoiseModel!=NULL) - { - if (dynamic_cast*>(m_NoiseModel)) - out.m_NoiseModel = new mitk::RicianNoiseModel(); - else if (dynamic_cast*>(m_NoiseModel)) - out.m_NoiseModel = new mitk::ChiSquareNoiseModel(); - out.m_NoiseModel->SetNoiseVariance(m_NoiseModel->GetNoiseVariance()); - } - - for (unsigned int i=0; i* outModel = NULL; - mitk::DiffusionSignalModel* signalModel = NULL; - 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 #include #include #include #include #include template< class ScalarType > mitk::FiberfoxParameters< ScalarType >::FiberfoxParameters() : m_NoiseModel(NULL) { } template< class ScalarType > mitk::FiberfoxParameters< ScalarType >::~FiberfoxParameters() { // if (m_NoiseModel!=NULL) // delete m_NoiseModel; } -void mitk::SignalGenerationParameters::GenerateGradientHalfShell() -{ - int NPoints = 2*m_NumGradients; - m_GradientDirections.clear(); - - m_NumBaseline = NPoints/20; - if (m_NumBaseline==0) - m_NumBaseline=1; - - GradientType g; - g.Fill(0.0); - for (unsigned int i=0; i theta; theta.set_size(NPoints); - vnl_vector phi; phi.set_size(NPoints); - double C = sqrt(4*M_PI); - phi(0) = 0.0; - phi(NPoints-1) = 0.0; - for(int i=0; i0 && i 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::SignalGenerationParameters::GradientType mitk::SignalGenerationParameters::GetGradientDirection(unsigned int i) -{ - return m_GradientDirections.at(i); -} - -void mitk::SignalGenerationParameters::SetNumWeightedVolumes(int numGradients) -{ - m_NumGradients = numGradients; - GenerateGradientHalfShell(); -} - -void mitk::SignalGenerationParameters::SetGradienDirections(GradientListType gradientList) -{ - m_GradientDirections = gradientList; - m_NumGradients = 0; - m_NumBaseline = 0; - for( unsigned int i=0; im_GradientDirections.size(); i++) - { - if (m_GradientDirections.at(i).GetNorm()>0.0001) - m_NumGradients++; - else - m_NumBaseline++; - } -} - -void mitk::SignalGenerationParameters::SetGradienDirections(mitk::DiffusionImage::GradientDirectionContainerType::Pointer gradientList) -{ - m_NumGradients = 0; - m_NumBaseline = 0; - m_GradientDirections.clear(); - - for( unsigned int i=0; iSize(); i++) - { - GradientType g; - g[0] = gradientList->at(i)[0]; - g[1] = gradientList->at(i)[1]; - g[2] = gradientList->at(i)[2]; - m_GradientDirections.push_back(g); - - if (m_GradientDirections.at(i).GetNorm()>0.0001) - m_NumGradients++; - else - m_NumBaseline++; - } -} - template< class ScalarType > -void mitk::FiberfoxParameters< ScalarType >::SaveParameters(string filename) +void mitk::FiberfoxParameters< ScalarType >::SaveParameters(std::string filename) { if(filename.empty()) return; if(".ffp"!=filename.substr(filename.size()-4, 4)) filename += ".ffp"; boost::property_tree::ptree parameters; // fiber generation parameters parameters.put("fiberfox.fibers.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.1", m_SignalGen.m_ImageDirection[0][0]); parameters.put("fiberfox.image.basic.direction.2", m_SignalGen.m_ImageDirection[0][1]); parameters.put("fiberfox.image.basic.direction.3", m_SignalGen.m_ImageDirection[0][2]); parameters.put("fiberfox.image.basic.direction.4", m_SignalGen.m_ImageDirection[1][0]); parameters.put("fiberfox.image.basic.direction.5", m_SignalGen.m_ImageDirection[1][1]); parameters.put("fiberfox.image.basic.direction.6", m_SignalGen.m_ImageDirection[1][2]); parameters.put("fiberfox.image.basic.direction.7", m_SignalGen.m_ImageDirection[2][0]); parameters.put("fiberfox.image.basic.direction.8", m_SignalGen.m_ImageDirection[2][1]); parameters.put("fiberfox.image.basic.direction.9", m_SignalGen.m_ImageDirection[2][2]); parameters.put("fiberfox.image.basic.numgradients", m_SignalGen.GetNumWeightedVolumes()); for( unsigned int i=0; im_SignalGen.GetNumVolumes(); i++) { - parameters.put("fiberfox.image.gradients."+boost::lexical_cast(i)+".x", m_SignalGen.GetGradientDirection(i)[0]); - parameters.put("fiberfox.image.gradients."+boost::lexical_cast(i)+".y", m_SignalGen.GetGradientDirection(i)[1]); - parameters.put("fiberfox.image.gradients."+boost::lexical_cast(i)+".z", m_SignalGen.GetGradientDirection(i)[2]); + parameters.put("fiberfox.image.gradients."+boost::lexical_cast(i)+".x", m_SignalGen.GetGradientDirection(i)[0]); + parameters.put("fiberfox.image.gradients."+boost::lexical_cast(i)+".y", m_SignalGen.GetGradientDirection(i)[1]); + parameters.put("fiberfox.image.gradients."+boost::lexical_cast(i)+".z", m_SignalGen.GetGradientDirection(i)[2]); } parameters.put("fiberfox.image.signalScale", m_SignalGen.m_SignalScale); parameters.put("fiberfox.image.tEcho", m_SignalGen.m_tEcho); parameters.put("fiberfox.image.tLine", m_SignalGen.m_tLine); parameters.put("fiberfox.image.tInhom", m_SignalGen.m_tInhom); parameters.put("fiberfox.image.bvalue", m_SignalGen.m_Bvalue); parameters.put("fiberfox.image.simulatekspace", m_SignalGen.m_SimulateKspaceAcquisition); parameters.put("fiberfox.image.axonRadius", m_SignalGen.m_AxonRadius); parameters.put("fiberfox.image.diffusiondirectionmode", m_SignalGen.m_DiffusionDirectionMode); parameters.put("fiberfox.image.fiberseparationthreshold", m_SignalGen.m_FiberSeparationThreshold); parameters.put("fiberfox.image.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.addringing", m_SignalGen.m_DoAddGibbsRinging); 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.addnoise", m_Misc.m_CheckAddNoiseBox); parameters.put("fiberfox.image.artifacts.addghosts", m_Misc.m_CheckAddGhostsBox); parameters.put("fiberfox.image.artifacts.addaliasing", m_Misc.m_CheckAddAliasingBox); parameters.put("fiberfox.image.artifacts.addspikes", m_Misc.m_CheckAddSpikesBox); parameters.put("fiberfox.image.artifacts.addeddycurrents", m_Misc.m_CheckAddEddyCurrentsBox); parameters.put("fiberfox.image.artifacts.doAddDistortions", m_Misc.m_CheckAddDistortionsBox); 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!=NULL) { parameters.put("fiberfox.image.artifacts.noisevariance", m_NoiseModel->GetNoiseVariance()); if (dynamic_cast*>(m_NoiseModel)) parameters.put("fiberfox.image.artifacts.noisetype", "rice"); else if (dynamic_cast*>(m_NoiseModel)) parameters.put("fiberfox.image.artifacts.noisetype", "chisquare"); } for (int i=0; i* signalModel = NULL; if (i(i)+".type", "fiber"); + parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".type", "fiber"); } else { signalModel = m_NonFiberModelList.at(i-m_FiberModelList.size()); - parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".type", "non-fiber"); + parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".type", "non-fiber"); } if (dynamic_cast*>(signalModel)) { mitk::StickModel* model = dynamic_cast*>(signalModel); - parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".stick.d", model->GetDiffusivity()); - parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".stick.t2", model->GetT2()); + parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".stick.d", model->GetDiffusivity()); + parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".stick.t2", model->GetT2()); } else if (dynamic_cast*>(signalModel)) { mitk::TensorModel* model = dynamic_cast*>(signalModel); - parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".tensor.d1", model->GetDiffusivity1()); - parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".tensor.d2", model->GetDiffusivity2()); - parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".tensor.d3", model->GetDiffusivity3()); - parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".tensor.t2", model->GetT2()); + parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".tensor.d1", model->GetDiffusivity1()); + parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".tensor.d2", model->GetDiffusivity2()); + parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".tensor.d3", model->GetDiffusivity3()); + parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".tensor.t2", model->GetT2()); } else if (dynamic_cast*>(signalModel)) { mitk::RawShModel* model = dynamic_cast*>(signalModel); - parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".prototype.minFA", model->GetFaRange().first); - parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".prototype.maxFA", model->GetFaRange().second); - parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".prototype.minADC", model->GetAdcRange().first); - parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".prototype.maxADC", model->GetAdcRange().second); - parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".prototype.numSamples", model->GetMaxNumKernels()); + parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".prototype.minFA", model->GetFaRange().first); + parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".prototype.maxFA", model->GetFaRange().second); + parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".prototype.minADC", model->GetAdcRange().first); + parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".prototype.maxADC", model->GetAdcRange().second); + parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".prototype.numSamples", model->GetMaxNumKernels()); } else if (dynamic_cast*>(signalModel)) { mitk::BallModel* model = dynamic_cast*>(signalModel); - parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".ball.d", model->GetDiffusivity()); - parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".ball.t2", model->GetT2()); + parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".ball.d", model->GetDiffusivity()); + parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".ball.t2", model->GetT2()); } else if (dynamic_cast*>(signalModel)) { mitk::AstroStickModel* model = dynamic_cast*>(signalModel); - parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".astrosticks.d", model->GetDiffusivity()); - parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".astrosticks.t2", model->GetT2()); - parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".astrosticks.randomize", model->GetRandomizeSticks()); + parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".astrosticks.d", model->GetDiffusivity()); + parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".astrosticks.t2", model->GetT2()); + parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".astrosticks.randomize", model->GetRandomizeSticks()); } else if (dynamic_cast*>(signalModel)) { mitk::DotModel* model = dynamic_cast*>(signalModel); - parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".dot.t2", model->GetT2()); + parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".dot.t2", model->GetT2()); } if (signalModel!=NULL) { - parameters.put("fiberfox.image.compartments."+boost::lexical_cast(i)+".ID", signalModel->m_CompartmentId); + parameters.put("fiberfox.image.compartments."+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)+".nrrd"); + writer->SetFileName(filename+"_VOLUME"+boost::lexical_cast(signalModel->m_CompartmentId)+".nrrd"); writer->SetInput(signalModel->GetVolumeFractionImage()); writer->Update(); - MITK_INFO << "Volume fraction image for compartment "+boost::lexical_cast(signalModel->m_CompartmentId)+" saved."; + 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.nrrd"); 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.nrrd"); writer->SetInput(m_SignalGen.m_MaskImage); writer->Update(); } catch(...) { MITK_INFO << "No mask image saved."; } } template< class ScalarType > -void mitk::FiberfoxParameters< ScalarType >::LoadParameters(string filename) +void mitk::FiberfoxParameters< ScalarType >::LoadParameters(std::string filename) { if(filename.empty()) return; boost::property_tree::ptree parameterTree; boost::property_tree::xml_parser::read_xml(filename, parameterTree); m_FiberModelList.clear(); m_NonFiberModelList.clear(); if (m_NoiseModel!=NULL) delete m_NoiseModel; BOOST_FOREACH( boost::property_tree::ptree::value_type const& v1, parameterTree.get_child("fiberfox") ) { if( v1.first == "fibers" ) { m_Misc.m_CheckRealTimeFibersBox = v1.second.get("realtime", m_Misc.m_CheckRealTimeFibersBox); m_Misc.m_CheckAdvancedFiberOptionsBox = v1.second.get("showadvanced", m_Misc.m_CheckAdvancedFiberOptionsBox); m_Misc.m_CheckConstantRadiusBox = v1.second.get("constantradius", m_Misc.m_CheckConstantRadiusBox); m_Misc.m_CheckIncludeFiducialsBox = v1.second.get("includeFiducials", m_Misc.m_CheckIncludeFiducialsBox); switch (v1.second.get("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 = v1.second.get("variance", m_FiberGen.m_Variance); m_FiberGen.m_Density = v1.second.get("density", m_FiberGen.m_Density); m_FiberGen.m_Sampling = v1.second.get("spline.sampling", m_FiberGen.m_Sampling); m_FiberGen.m_Tension = v1.second.get("spline.tension", m_FiberGen.m_Tension); m_FiberGen.m_Continuity = v1.second.get("spline.continuity", m_FiberGen.m_Continuity); m_FiberGen.m_Bias = v1.second.get("spline.bias", m_FiberGen.m_Bias); m_FiberGen.m_Rotation[0] = v1.second.get("rotation.x", m_FiberGen.m_Rotation[0]); m_FiberGen.m_Rotation[1] = v1.second.get("rotation.y", m_FiberGen.m_Rotation[1]); m_FiberGen.m_Rotation[2] = v1.second.get("rotation.z", m_FiberGen.m_Rotation[2]); m_FiberGen.m_Translation[0] = v1.second.get("translation.x", m_FiberGen.m_Translation[0]); m_FiberGen.m_Translation[1] = v1.second.get("translation.y", m_FiberGen.m_Translation[1]); m_FiberGen.m_Translation[2] = v1.second.get("translation.z", m_FiberGen.m_Translation[2]); m_FiberGen.m_Scale[0] = v1.second.get("scale.x", m_FiberGen.m_Scale[0]); m_FiberGen.m_Scale[1] = v1.second.get("scale.y", m_FiberGen.m_Scale[1]); m_FiberGen.m_Scale[2] = v1.second.get("scale.z", m_FiberGen.m_Scale[2]); } else if ( v1.first == "image" ) { - m_Misc.m_SignalModelString = v1.second.get("signalmodelstring", m_Misc.m_SignalModelString); - m_Misc.m_ArtifactModelString = v1.second.get("artifactmodelstring", m_Misc.m_ArtifactModelString); - m_Misc.m_OutputPath = v1.second.get("outpath", m_Misc.m_OutputPath); + m_Misc.m_SignalModelString = v1.second.get("signalmodelstring", m_Misc.m_SignalModelString); + m_Misc.m_ArtifactModelString = v1.second.get("artifactmodelstring", m_Misc.m_ArtifactModelString); + m_Misc.m_OutputPath = v1.second.get("outpath", m_Misc.m_OutputPath); m_Misc.m_CheckOutputVolumeFractionsBox = v1.second.get("outputvolumefractions", m_Misc.m_CheckOutputVolumeFractionsBox); m_Misc.m_CheckAdvancedSignalOptionsBox = v1.second.get("showadvanced", m_Misc.m_CheckAdvancedSignalOptionsBox); m_Misc.m_CheckAddDistortionsBox = v1.second.get("artifacts.doAddDistortions", m_Misc.m_CheckAddDistortionsBox); m_Misc.m_CheckAddNoiseBox = v1.second.get("artifacts.addnoise", m_Misc.m_CheckAddNoiseBox); m_Misc.m_CheckAddGhostsBox = v1.second.get("artifacts.addghosts", m_Misc.m_CheckAddGhostsBox); m_Misc.m_CheckAddAliasingBox = v1.second.get("artifacts.addaliasing", m_Misc.m_CheckAddAliasingBox); m_Misc.m_CheckAddSpikesBox = v1.second.get("artifacts.addspikes", m_Misc.m_CheckAddSpikesBox); m_Misc.m_CheckAddEddyCurrentsBox = v1.second.get("artifacts.addeddycurrents", m_Misc.m_CheckAddEddyCurrentsBox); m_SignalGen.m_ImageRegion.SetSize(0, v1.second.get("basic.size.x",m_SignalGen.m_ImageRegion.GetSize(0))); m_SignalGen.m_ImageRegion.SetSize(1, v1.second.get("basic.size.y",m_SignalGen.m_ImageRegion.GetSize(1))); m_SignalGen.m_ImageRegion.SetSize(2, v1.second.get("basic.size.z",m_SignalGen.m_ImageRegion.GetSize(2))); m_SignalGen.m_ImageSpacing[0] = v1.second.get("basic.spacing.x",m_SignalGen.m_ImageSpacing[0]); m_SignalGen.m_ImageSpacing[1] = v1.second.get("basic.spacing.y",m_SignalGen.m_ImageSpacing[1]); m_SignalGen.m_ImageSpacing[2] = v1.second.get("basic.spacing.z",m_SignalGen.m_ImageSpacing[2]); m_SignalGen.m_ImageOrigin[0] = v1.second.get("basic.origin.x",m_SignalGen.m_ImageOrigin[0]); m_SignalGen.m_ImageOrigin[1] = v1.second.get("basic.origin.y",m_SignalGen.m_ImageOrigin[1]); m_SignalGen.m_ImageOrigin[2] = v1.second.get("basic.origin.z",m_SignalGen.m_ImageOrigin[2]); m_SignalGen.m_ImageDirection[0][0] = v1.second.get("basic.direction.1",m_SignalGen.m_ImageDirection[0][0]); m_SignalGen.m_ImageDirection[0][1] = v1.second.get("basic.direction.2",m_SignalGen.m_ImageDirection[0][1]); m_SignalGen.m_ImageDirection[0][2] = v1.second.get("basic.direction.3",m_SignalGen.m_ImageDirection[0][2]); m_SignalGen.m_ImageDirection[1][0] = v1.second.get("basic.direction.4",m_SignalGen.m_ImageDirection[1][0]); m_SignalGen.m_ImageDirection[1][1] = v1.second.get("basic.direction.5",m_SignalGen.m_ImageDirection[1][1]); m_SignalGen.m_ImageDirection[1][2] = v1.second.get("basic.direction.6",m_SignalGen.m_ImageDirection[1][2]); m_SignalGen.m_ImageDirection[2][0] = v1.second.get("basic.direction.7",m_SignalGen.m_ImageDirection[2][0]); m_SignalGen.m_ImageDirection[2][1] = v1.second.get("basic.direction.8",m_SignalGen.m_ImageDirection[2][1]); m_SignalGen.m_ImageDirection[2][2] = v1.second.get("basic.direction.9",m_SignalGen.m_ImageDirection[2][2]); m_SignalGen.m_SignalScale = v1.second.get("signalScale", m_SignalGen.m_SignalScale); m_SignalGen.m_tEcho = v1.second.get("tEcho", m_SignalGen.m_tEcho); m_SignalGen.m_tLine = v1.second.get("tLine", m_SignalGen.m_tLine); m_SignalGen.m_tInhom = v1.second.get("tInhom", m_SignalGen.m_tInhom); m_SignalGen.m_Bvalue = v1.second.get("bvalue", m_SignalGen.m_Bvalue); m_SignalGen.m_SimulateKspaceAcquisition = v1.second.get("simulatekspace", m_SignalGen.m_SimulateKspaceAcquisition); m_SignalGen.m_AxonRadius = v1.second.get("axonRadius", m_SignalGen.m_AxonRadius); switch (v1.second.get("diffusiondirectionmode", 0)) { case 0: m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::FIBER_TANGENT_DIRECTIONS; break; case 1: m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::MAIN_FIBER_DIRECTIONS; break; case 2: m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::RANDOM_DIRECTIONS; break; default: m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::FIBER_TANGENT_DIRECTIONS; } m_SignalGen.m_FiberSeparationThreshold = v1.second.get("fiberseparationthreshold", m_SignalGen.m_FiberSeparationThreshold); m_SignalGen.m_Spikes = v1.second.get("artifacts.spikesnum", m_SignalGen.m_Spikes); m_SignalGen.m_SpikeAmplitude = v1.second.get("artifacts.spikesscale", m_SignalGen.m_SpikeAmplitude); m_SignalGen.m_KspaceLineOffset = v1.second.get("artifacts.kspaceLineOffset", m_SignalGen.m_KspaceLineOffset); m_SignalGen.m_EddyStrength = v1.second.get("artifacts.eddyStrength", m_SignalGen.m_EddyStrength); m_SignalGen.m_Tau = v1.second.get("artifacts.eddyTau", m_SignalGen.m_Tau); m_SignalGen.m_CroppingFactor = v1.second.get("artifacts.aliasingfactor", m_SignalGen.m_CroppingFactor); m_SignalGen.m_DoAddGibbsRinging = v1.second.get("artifacts.addringing", m_SignalGen.m_DoAddGibbsRinging); m_SignalGen.m_DoSimulateRelaxation = v1.second.get("doSimulateRelaxation", m_SignalGen.m_DoSimulateRelaxation); m_SignalGen.m_DoDisablePartialVolume = v1.second.get("doDisablePartialVolume", m_SignalGen.m_DoDisablePartialVolume); m_SignalGen.m_DoAddMotion = v1.second.get("artifacts.doAddMotion", m_SignalGen.m_DoAddMotion); m_SignalGen.m_DoRandomizeMotion = v1.second.get("artifacts.randomMotion", m_SignalGen.m_DoRandomizeMotion); m_SignalGen.m_Translation[0] = v1.second.get("artifacts.translation0", m_SignalGen.m_Translation[0]); m_SignalGen.m_Translation[1] = v1.second.get("artifacts.translation1", m_SignalGen.m_Translation[1]); m_SignalGen.m_Translation[2] = v1.second.get("artifacts.translation2", m_SignalGen.m_Translation[2]); m_SignalGen.m_Rotation[0] = v1.second.get("artifacts.rotation0", m_SignalGen.m_Rotation[0]); m_SignalGen.m_Rotation[1] = v1.second.get("artifacts.rotation1", m_SignalGen.m_Rotation[1]); m_SignalGen.m_Rotation[2] = v1.second.get("artifacts.rotation2", m_SignalGen.m_Rotation[2]); // m_SignalGen.SetNumWeightedVolumes(v1.second.get("numgradients", m_SignalGen.GetNumWeightedVolumes())); SignalGenerationParameters::GradientListType gradients; BOOST_FOREACH( boost::property_tree::ptree::value_type const& v2, v1.second.get_child("gradients") ) { SignalGenerationParameters::GradientType g; g[0] = v2.second.get("x"); g[1] = v2.second.get("y"); g[2] = v2.second.get("z"); gradients.push_back(g); } m_SignalGen.SetGradienDirections(gradients); try { - if (v1.second.get("artifacts.noisetype")=="rice") + if (v1.second.get("artifacts.noisetype")=="rice") { m_NoiseModel = new mitk::RicianNoiseModel(); m_NoiseModel->SetNoiseVariance(v1.second.get("artifacts.noisevariance")); } } catch(...) {} try { - if (v1.second.get("artifacts.noisetype")=="chisquare") + if (v1.second.get("artifacts.noisetype")=="chisquare") { m_NoiseModel = new mitk::ChiSquareNoiseModel(); m_NoiseModel->SetNoiseVariance(v1.second.get("artifacts.noisevariance")); } } catch(...){ } BOOST_FOREACH( boost::property_tree::ptree::value_type const& v2, v1.second.get_child("compartments") ) { mitk::DiffusionSignalModel* signalModel = NULL; try // stick model { mitk::StickModel* model = new mitk::StickModel(); model->SetDiffusivity(v2.second.get("stick.d")); model->SetT2(v2.second.get("stick.t2")); model->m_CompartmentId = v2.second.get("ID"); - if (v2.second.get("type")=="fiber") + if (v2.second.get("type")=="fiber") m_FiberModelList.push_back(model); - else if (v2.second.get("type")=="non-fiber") + else if (v2.second.get("type")=="non-fiber") m_NonFiberModelList.push_back(model); signalModel = model; } catch (...) { } try // tensor model { mitk::TensorModel* model = new mitk::TensorModel(); model->SetDiffusivity1(v2.second.get("tensor.d1")); model->SetDiffusivity2(v2.second.get("tensor.d2")); model->SetDiffusivity3(v2.second.get("tensor.d3")); model->SetT2(v2.second.get("tensor.t2")); model->m_CompartmentId = v2.second.get("ID"); - if (v2.second.get("type")=="fiber") + if (v2.second.get("type")=="fiber") m_FiberModelList.push_back(model); - else if (v2.second.get("type")=="non-fiber") + else if (v2.second.get("type")=="non-fiber") m_NonFiberModelList.push_back(model); signalModel = model; } catch (...) { } try // ball model { mitk::BallModel* model = new mitk::BallModel(); model->SetDiffusivity(v2.second.get("ball.d")); model->SetT2(v2.second.get("ball.t2")); model->m_CompartmentId = v2.second.get("ID"); - if (v2.second.get("type")=="fiber") + if (v2.second.get("type")=="fiber") m_FiberModelList.push_back(model); - else if (v2.second.get("type")=="non-fiber") + else if (v2.second.get("type")=="non-fiber") m_NonFiberModelList.push_back(model); signalModel = model; } catch (...) { } try // astrosticks model { mitk::AstroStickModel* model = new AstroStickModel(); model->SetDiffusivity(v2.second.get("astrosticks.d")); model->SetT2(v2.second.get("astrosticks.t2")); model->SetRandomizeSticks(v2.second.get("astrosticks.randomize")); model->m_CompartmentId = v2.second.get("ID"); - if (v2.second.get("type")=="fiber") + if (v2.second.get("type")=="fiber") m_FiberModelList.push_back(model); - else if (v2.second.get("type")=="non-fiber") + else if (v2.second.get("type")=="non-fiber") m_NonFiberModelList.push_back(model); signalModel = model; } catch (...) { } try // dot model { mitk::DotModel* model = new mitk::DotModel(); model->SetT2(v2.second.get("dot.t2")); model->m_CompartmentId = v2.second.get("ID"); - if (v2.second.get("type")=="fiber") + if (v2.second.get("type")=="fiber") m_FiberModelList.push_back(model); - else if (v2.second.get("type")=="non-fiber") + else if (v2.second.get("type")=="non-fiber") m_NonFiberModelList.push_back(model); signalModel = model; } catch (...) { } try // prototype model { mitk::RawShModel* model = new mitk::RawShModel(); model->SetMaxNumKernels(v2.second.get("prototype.numSamples")); model->SetFaRange(v2.second.get("prototype.minFA"), v2.second.get("prototype.maxFA")); model->SetAdcRange(v2.second.get("prototype.minADC"), v2.second.get("prototype.maxADC")); model->m_CompartmentId = v2.second.get("ID"); - if (v2.second.get("type")=="fiber") + if (v2.second.get("type")=="fiber") m_FiberModelList.push_back(model); - else if (v2.second.get("type")=="non-fiber") + else if (v2.second.get("type")=="non-fiber") m_NonFiberModelList.push_back(model); signalModel = model; } catch (...) { } if (signalModel!=NULL) { signalModel->SetGradientList(gradients); try{ itk::ImageFileReader::Pointer reader = itk::ImageFileReader::New(); - reader->SetFileName(filename+"_VOLUME"+v2.second.get("ID")+".nrrd"); + reader->SetFileName(filename+"_VOLUME"+v2.second.get("ID")+".nrrd"); reader->Update(); signalModel->SetVolumeFractionImage(reader->GetOutput()); } catch(...) { } } } } } try{ itk::ImageFileReader::Pointer reader = itk::ImageFileReader::New(); reader->SetFileName(filename+"_FMAP.nrrd"); reader->Update(); m_SignalGen.m_FrequencyMap = reader->GetOutput(); } catch(...) { MITK_INFO << "No frequency map saved."; } try{ itk::ImageFileReader::Pointer reader = itk::ImageFileReader::New(); reader->SetFileName(filename+"_MASK.nrrd"); reader->Update(); m_SignalGen.m_MaskImage = reader->GetOutput(); } catch(...) { MITK_INFO << "No mask image saved."; } } template< class ScalarType > void mitk::FiberfoxParameters< ScalarType >::PrintSelf() { MITK_INFO << "Not implemented :("; } diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h index 4e43573b4b..7d110c4b8f 100644 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberfoxParameters.h @@ -1,321 +1,122 @@ /*=================================================================== 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 -using namespace std; - namespace mitk { -/** Signal generation */ -class MitkFiberTracking_EXPORT SignalGenerationParameters -{ -public: - typedef itk::Image ItkDoubleImgType; - typedef itk::Image ItkUcharImgType; - typedef itk::Vector GradientType; - typedef std::vector GradientListType; - - enum DiffusionDirectionMode { - FIBER_TANGENT_DIRECTIONS, - MAIN_FIBER_DIRECTIONS, - RANDOM_DIRECTIONS - }; - - SignalGenerationParameters() - : m_FrequencyMap(NULL) - , m_MaskImage(NULL) - , m_SignalScale(100) - , m_tEcho(100) - , m_tLine(1) - , m_tInhom(50) - , m_Bvalue(1000) - , m_SimulateKspaceAcquisition(false) - , m_AxonRadius(0) - , m_DiffusionDirectionMode(SignalGenerationParameters::FIBER_TANGENT_DIRECTIONS) - , m_FiberSeparationThreshold(30) - , m_Spikes(0) - , m_SpikeAmplitude(1) - , m_KspaceLineOffset(0) - , m_EddyStrength(0) - , m_Tau(70) - , m_CroppingFactor(1) - , m_DoAddGibbsRinging(false) - , m_DoSimulateRelaxation(true) - , m_DoDisablePartialVolume(false) - , m_DoAddMotion(false) - , m_DoRandomizeMotion(true) - { - 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_ImageRegion; ///< Image size. - itk::Vector m_ImageSpacing; ///< Image voxel size. - itk::Point m_ImageOrigin; ///< Image origin. - itk::Matrix m_ImageDirection; ///< Image rotation matrix. - - /** Other acquisitions parameters */ - double m_SignalScale; ///< Scaling factor for output signal (before noise is added). - double m_tEcho; ///< Echo time TE. - double m_tLine; ///< k-space line readout time. - double m_tInhom; ///< T2' - double m_Bvalue; ///< Acquisition b-value - bool m_SimulateKspaceAcquisition;///< 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. - DiffusionDirectionMode m_DiffusionDirectionMode; ///< Determines how the main diffusion direction of the signal models is selected - double m_FiberSeparationThreshold; ///< Used for random and and mein fiber deriction DiffusionDirectionMode - - /** Artifacts and other effects */ - unsigned int m_Spikes; ///< Number of spikes randomly appearing in the image - double m_SpikeAmplitude; ///< amplitude of spikes relative to the largest signal intensity (magnitude of complex) - double m_KspaceLineOffset; ///< Causes N/2 ghosts. Larger offset means stronger ghost. - double m_EddyStrength; ///< Strength of eddy current induced gradients in mT/m. - double m_Tau; ///< Eddy current decay constant (in ms) - double m_CroppingFactor; ///< FOV size in y-direction is multiplied by this factor. Causes aliasing artifacts. - 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. - itk::Vector m_Translation; ///< Maximum translational motion. - itk::Vector m_Rotation; ///< Maximum rotational motion. - ItkDoubleImgType::Pointer m_FrequencyMap; ///< If != NULL, distortions are added to the image using this frequency map. - ItkUcharImgType::Pointer m_MaskImage; ///< Signal is only genrated inside of the mask image. - - inline void GenerateGradientHalfShell(); ///< Generates half shell of gradient directions (with m_NumGradients non-zero directions) - inline std::vector< int > GetBaselineIndices(); ///< Returns list of nun-diffusion-weighted image volume indices - inline unsigned int GetFirstBaselineIndex(); ///< Returns index of first non-diffusion-weighted image volume - inline bool IsBaselineIndex(unsigned int idx); ///< Checks if image volume with given index is non-diffusion-weighted volume or not. - inline unsigned int GetNumWeightedVolumes(); ///< Get number of diffusion-weighted image volumes - inline unsigned int GetNumBaselineVolumes(); ///< Get number of non-diffusion-weighted image volumes - inline unsigned int GetNumVolumes(); ///< Get number of baseline and diffusion-weighted image volumes - inline GradientListType GetGradientDirections(); ///< Return gradient direction container - inline GradientType GetGradientDirection(unsigned int i); - - inline void SetNumWeightedVolumes(int numGradients); ///< Automaticall calls GenerateGradientHalfShell() afterwards. - inline void SetGradienDirections(GradientListType gradientList); - inline void SetGradienDirections(mitk::DiffusionImage::GradientDirectionContainerType::Pointer gradientList); - -protected: - - unsigned int m_NumGradients; ///< Number of diffusion-weighted image volumes. - unsigned int m_NumBaseline; ///< Number of non-diffusion-weighted image volumes. - GradientListType m_GradientDirections; ///< Total number of image volumes. -}; - -/** 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 vector< vector< mitk::PlanarEllipse::Pointer > > FiducialListType; - typedef vector< 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(NULL) - , m_SignalModelString("") - , m_ArtifactModelString("") - , m_OutputPath("") - , m_CheckOutputVolumeFractionsBox(false) - , m_CheckAdvancedSignalOptionsBox(false) - , m_CheckAddNoiseBox(false) - , m_CheckAddGhostsBox(false) - , m_CheckAddAliasingBox(false) - , m_CheckAddSpikesBox(false) - , m_CheckAddEddyCurrentsBox(false) - , m_CheckAddDistortionsBox(false) - , 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. - string m_SignalModelString; ///< Appendet to the name of the result node - string m_ArtifactModelString; ///< Appendet to the name of the result node - string m_OutputPath; ///< Image is automatically saved to the specified folder after simulation is finished. - - /** member variables that store the check-state of GUI checkboxes */ - // image generation - bool m_CheckOutputVolumeFractionsBox; - bool m_CheckAdvancedSignalOptionsBox; - bool m_CheckAddNoiseBox; - bool m_CheckAddGhostsBox; - bool m_CheckAddAliasingBox; - bool m_CheckAddSpikesBox; - bool m_CheckAddEddyCurrentsBox; - bool m_CheckAddDistortionsBox; - // fiber generation - bool m_CheckRealTimeFibersBox; - bool m_CheckAdvancedFiberOptionsBox; - bool m_CheckConstantRadiusBox; - bool m_CheckIncludeFiducialsBox; -}; - /** * \brief Datastructure to manage the Fiberfox signal generation parameters. * */ template< class ScalarType = double > -class MitkFiberTracking_EXPORT FiberfoxParameters +class FiberfoxParameters { public: typedef itk::Image ItkDoubleImgType; typedef itk::Image ItkUcharImgType; typedef DiffusionSignalModel DiffusionModelType; typedef std::vector< DiffusionModelType* > DiffusionModelListType; typedef DiffusionNoiseModel NoiseModelType; FiberfoxParameters(); ~FiberfoxParameters(); /** Get same parameter object with different template parameter */ template< class OutType > FiberfoxParameters< OutType > CopyParameters() { FiberfoxParameters< OutType > out; out.m_FiberGen = m_FiberGen; out.m_SignalGen = m_SignalGen; out.m_Misc = m_Misc; if (m_NoiseModel!=NULL) { if (dynamic_cast*>(m_NoiseModel)) out.m_NoiseModel = new mitk::RicianNoiseModel(); else if (dynamic_cast*>(m_NoiseModel)) out.m_NoiseModel = new mitk::ChiSquareNoiseModel(); out.m_NoiseModel->SetNoiseVariance(m_NoiseModel->GetNoiseVariance()); } for (unsigned int i=0; i* outModel = NULL; mitk::DiffusionSignalModel* signalModel = NULL; 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 #include #include #include #include #include +#include +#include +#include using namespace mitk; using namespace boost::math; template< class ScalarType > RawShModel< ScalarType >::RawShModel() : m_ShOrder(0) , m_ModelIndex(-1) , m_MaxNumKernels(1000) { this->m_RandGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); this->m_RandGen->SetSeed(); m_AdcRange.first = 0; m_AdcRange.second = 0.004; m_FaRange.first = 0; m_FaRange.second = 1; } template< class ScalarType > RawShModel< ScalarType >::~RawShModel() { } template< class ScalarType > void RawShModel< ScalarType >::Clear() { m_ShCoefficients.clear(); m_PrototypeMaxDirection.clear(); m_B0Signal.clear(); } template< class ScalarType > void RawShModel< ScalarType >::RandomModel() { m_ModelIndex = this->m_RandGen->GetIntegerVariate(m_B0Signal.size()-1); } template< class ScalarType > unsigned int RawShModel< ScalarType >::GetNumberOfKernels() { return m_B0Signal.size(); } +template< class ScalarType > +bool RawShModel< ScalarType >::SampleKernels(DiffusionImage::Pointer diffImg, ItkUcharImageType::Pointer maskImage, TensorImageType::Pointer tensorImage, QballFilterType::CoefficientImageType::Pointer itkFeatureImage, ItkDoubleImageType::Pointer adcImage) +{ + if (diffImg.IsNull()) + return false; + + const int shOrder = 2; + + if (tensorImage.IsNull()) + { + typedef itk::DiffusionTensor3DReconstructionImageFilter< short, short, double > TensorReconstructionImageFilterType; + TensorReconstructionImageFilterType::Pointer filter = TensorReconstructionImageFilterType::New(); + filter->SetGradientImage( diffImg->GetDirections(), diffImg->GetVectorImage() ); + filter->SetBValue(diffImg->GetReferenceBValue()); + filter->Update(); + tensorImage = filter->GetOutput(); + } + + const int NumCoeffs = (shOrder*shOrder + shOrder + 2)/2 + shOrder; + if (itkFeatureImage.IsNull()) + { + QballFilterType::Pointer qballfilter = QballFilterType::New(); + qballfilter->SetGradientImage( diffImg->GetDirections(), diffImg->GetVectorImage() ); + qballfilter->SetBValue(diffImg->GetReferenceBValue()); + qballfilter->SetLambda(0.006); + qballfilter->SetNormalizationMethod(QballFilterType::QBAR_RAW_SIGNAL); + qballfilter->Update(); + itkFeatureImage = qballfilter->GetCoefficientImage(); + } + + if (adcImage.IsNull()) + { + itk::AdcImageFilter< short, double >::Pointer adcFilter = itk::AdcImageFilter< short, double >::New(); + adcFilter->SetInput(diffImg->GetVectorImage()); + adcFilter->SetGradientDirections(diffImg->GetDirections()); + adcFilter->SetB_value(diffImg->GetReferenceBValue()); + adcFilter->Update(); + adcImage = adcFilter->GetOutput(); + } + + int b0Index; + for (unsigned int i=0; iGetDirectionsWithoutMeasurementFrame()->Size(); i++) + if ( diffImg->GetDirectionsWithoutMeasurementFrame()->GetElement(i).magnitude()<0.001 ) + { + b0Index = i; + break; + } + + double max = 0; + { + itk::ImageRegionIterator< itk::VectorImage< short, 3 > > it(diffImg->GetVectorImage(), diffImg->GetVectorImage()->GetLargestPossibleRegion()); + while(!it.IsAtEnd()) + { + if (maskImage.IsNotNull() && maskImage->GetPixel(it.GetIndex())<=0) + { + ++it; + continue; + } + if (it.Get()[b0Index]>max) + max = it.Get()[b0Index]; + ++it; + } + } + + MITK_INFO << "Sampling signal kernels."; + unsigned int count = 0; + itk::ImageRegionIterator< itk::Image< itk::DiffusionTensor3D< double >, 3 > > it(tensorImage, tensorImage->GetLargestPossibleRegion()); + while(!it.IsAtEnd()) + { + bool skipPixel = false; + if (maskImage.IsNotNull() && maskImage->GetPixel(it.GetIndex())<=0) + { + ++it; + continue; + } + for (unsigned int i=0; iGetDirections()->Size(); i++) + { + if (diffImg->GetVectorImage()->GetPixel(it.GetIndex())[i]!=diffImg->GetVectorImage()->GetPixel(it.GetIndex())[i] || diffImg->GetVectorImage()->GetPixel(it.GetIndex())[i]<=0 || diffImg->GetVectorImage()->GetPixel(it.GetIndex())[i]>diffImg->GetVectorImage()->GetPixel(it.GetIndex())[b0Index]) + { + skipPixel = true; + break; + } + } + if (skipPixel) + { + ++it; + continue; + } + + typedef itk::DiffusionTensor3D TensorType; + TensorType::EigenValuesArrayType eigenvalues; + TensorType::EigenVectorsMatrixType eigenvectors; + TensorType tensor = it.Get(); + double FA = tensor.GetFractionalAnisotropy(); + double ADC = adcImage->GetPixel(it.GetIndex()); + QballFilterType::CoefficientImageType::PixelType itkv = itkFeatureImage->GetPixel(it.GetIndex()); + vnl_vector_fixed< double, NumCoeffs > coeffs; + for (unsigned int c=0; cGetMaxNumKernels()>this->GetNumberOfKernels() && FA>m_FaRange.first && FAm_AdcRange.first && ADCSetShCoefficients( coeffs, (double)diffImg->GetVectorImage()->GetPixel(it.GetIndex())[b0Index]/max )) + { + tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); + itk::Vector dir; + dir[0] = eigenvectors(2, 0); + dir[1] = eigenvectors(2, 1); + dir[2] = eigenvectors(2, 2); + m_PrototypeMaxDirection.push_back(dir); + count++; + MITK_INFO << "KERNEL: " << it.GetIndex() << " (" << count << ")"; + } + } + ++it; + } + + if (m_ShCoefficients.size()>0) + return true; + return false; +} + // convert cartesian to spherical coordinates template< class ScalarType > void RawShModel< ScalarType >::Cart2Sph( GradientListType gradients ) { m_SphCoords.set_size(gradients.size(), 2); m_SphCoords.fill(0.0); for (unsigned int i=0; i 0.0001 ) { gradients[i].Normalize(); m_SphCoords(i,0) = acos(gradients[i][2]); // theta m_SphCoords(i,1) = atan2(gradients[i][1], gradients[i][0]); // phi } } } template< class ScalarType > void RawShModel< ScalarType >::SetFiberDirection(GradientType fiberDirection) { this->m_FiberDirection = fiberDirection; this->m_FiberDirection.Normalize(); RandomModel(); GradientType axis = itk::CrossProduct(this->m_FiberDirection, m_PrototypeMaxDirection.at(m_ModelIndex)); axis.Normalize(); vnl_quaternion rotation(axis.GetVnlVector(), acos(dot_product(this->m_FiberDirection.GetVnlVector(), m_PrototypeMaxDirection.at(m_ModelIndex).GetVnlVector()))); rotation.normalize(); GradientListType gradients; for (unsigned int i=0; im_GradientList.size(); i++) { GradientType dir = this->m_GradientList.at(i); if( dir.GetNorm() > 0.0001 ) { dir.Normalize(); vnl_vector_fixed< double, 3 > vnlDir = rotation.rotate(dir.GetVnlVector()); dir[0] = vnlDir[0]; dir[1] = vnlDir[1]; dir[2] = vnlDir[2]; dir.Normalize(); } gradients.push_back(dir); } Cart2Sph( gradients ); } template< class ScalarType > bool RawShModel< ScalarType >::SetShCoefficients(vnl_vector< double > shCoefficients, double b0 ) { m_ShOrder = 2; while ( (m_ShOrder*m_ShOrder + m_ShOrder + 2)/2 + m_ShOrder <= shCoefficients.size() ) m_ShOrder += 2; m_ShOrder -= 2; m_ModelIndex = m_B0Signal.size(); m_B0Signal.push_back(b0); m_ShCoefficients.push_back(shCoefficients); -// itk::OrientationDistributionFunction odf; -// GradientListType gradients; -// for (unsigned int i=0; im_GradientList ); -// PixelType signal = SimulateMeasurement(); - -// int minDirIdx = 0; -// double min = itk::NumericTraits::max(); -// for (unsigned int i=0; ib0 || signal[i]<0) -// { -// MITK_INFO << "Corrupted signal value detected. Kernel rejected."; -// m_B0Signal.pop_back(); -// m_ShCoefficients.pop_back(); -// return false; -// } -// if (signal[i]m_GradientList.at(minDirIdx); -// maxDir.Normalize(); -// m_PrototypeMaxDirection.push_back(maxDir); + // itk::OrientationDistributionFunction odf; + // GradientListType gradients; + // for (unsigned int i=0; im_GradientList ); + // PixelType signal = SimulateMeasurement(); + + // int minDirIdx = 0; + // double min = itk::NumericTraits::max(); + // for (unsigned int i=0; ib0 || signal[i]<0) + // { + // MITK_INFO << "Corrupted signal value detected. Kernel rejected."; + // m_B0Signal.pop_back(); + // m_ShCoefficients.pop_back(); + // return false; + // } + // if (signal[i]m_GradientList.at(minDirIdx); + // maxDir.Normalize(); + // m_PrototypeMaxDirection.push_back(maxDir); Cart2Sph( this->m_GradientList ); m_ModelIndex = -1; return true; } template< class ScalarType > ScalarType RawShModel< ScalarType >::SimulateMeasurement(unsigned int dir) { ScalarType signal = 0; if (m_ModelIndex==-1) RandomModel(); if (dir>=this->m_GradientList.size()) return signal; int j, m; double mag, plm; if (this->m_GradientList[dir].GetNorm()>0.001) { j=0; for (int l=0; l<=m_ShOrder; l=l+2) for (m=-l; m<=l; m++) { plm = legendre_p(l,abs(m),cos(m_SphCoords(dir,0))); mag = sqrt((double)(2*l+1)/(4.0*M_PI)*factorial(l-abs(m))/factorial(l+abs(m)))*plm; double basis; if (m<0) basis = sqrt(2.0)*mag*cos(fabs((double)m)*m_SphCoords(dir,1)); else if (m==0) basis = mag; else basis = pow(-1.0, m)*sqrt(2.0)*mag*sin(m*m_SphCoords(dir,1)); signal += m_ShCoefficients.at(m_ModelIndex)[j]*basis; j++; } } else signal = m_B0Signal.at(m_ModelIndex); m_ModelIndex = -1; return signal; } template< class ScalarType > typename RawShModel< ScalarType >::PixelType RawShModel< ScalarType >::SimulateMeasurement() { if (m_ModelIndex==-1) RandomModel(); PixelType signal; signal.SetSize(this->m_GradientList.size()); int M = this->m_GradientList.size(); int j, m; double mag, plm; vnl_matrix< double > shBasis; shBasis.set_size(M, m_ShCoefficients.at(m_ModelIndex).size()); shBasis.fill(0.0); for (int p=0; pm_GradientList[p].GetNorm()>0.001) { j=0; for (int l=0; l<=m_ShOrder; l=l+2) for (m=-l; m<=l; m++) { plm = legendre_p(l,abs(m),cos(m_SphCoords(p,0))); mag = sqrt((double)(2*l+1)/(4.0*M_PI)*factorial(l-abs(m))/factorial(l+abs(m)))*plm; if (m<0) shBasis(p,j) = sqrt(2.0)*mag*cos(fabs((double)m)*m_SphCoords(p,1)); else if (m==0) shBasis(p,j) = mag; else shBasis(p,j) = pow(-1.0, m)*sqrt(2.0)*mag*sin(m*m_SphCoords(p,1)); j++; } double val = 0; for (unsigned int cidx=0; cidx +#include +#include namespace mitk { /** * \brief The spherical harmonic representation of a prototype diffusion weighted MR signal is used to obtain the direction dependent signal. * */ template< class ScalarType = double > class RawShModel : public DiffusionSignalModel< ScalarType > { public: RawShModel(); template< class OtherType >RawShModel(RawShModel* model) { this->m_CompartmentId = model->m_CompartmentId; this->m_T2 = model->GetT2(); this->m_FiberDirection = model->GetFiberDirection(); this->m_GradientList = model->GetGradientList(); this->m_VolumeFractionImage = model->GetVolumeFractionImage(); this->m_RandGen = model->GetRandomGenerator(); this->m_AdcRange = model->GetAdcRange(); this->m_FaRange = model->GetFaRange(); this->m_ShCoefficients = model->GetShCoefficients(); this->m_B0Signal = model->GetB0Signals(); this->m_SphCoords = model->GetSphericalCoordinates(); this->m_ShOrder = model->GetShOrder(); this->m_ModelIndex = model->GetModelIndex(); this->m_MaxNumKernels = model->GetMaxNumKernels(); } ~RawShModel(); + typedef itk::Image< double, 3 > ItkDoubleImageType; + typedef itk::Image< unsigned char, 3 > ItkUcharImageType; + typedef itk::Image< itk::DiffusionTensor3D< double >, 3 > TensorImageType; + typedef itk::AnalyticalDiffusionQballReconstructionImageFilter QballFilterType; typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; typedef itk::Matrix< double, 3, 3 > MatrixType; /** Actual signal generation **/ PixelType SimulateMeasurement(); ScalarType SimulateMeasurement(unsigned int dir); bool SetShCoefficients(vnl_vector< double > shCoefficients, double b0); void SetFiberDirection(GradientType fiberDirection); void SetGradientList(GradientListType gradientList) { this->m_GradientList = gradientList; } void SetFaRange(double min, double max){ m_FaRange.first = min; m_FaRange.second = max; } void SetAdcRange(double min, double max){ m_AdcRange.first = min; m_AdcRange.second = max; } void SetMaxNumKernels(unsigned int max){ m_MaxNumKernels = max; } unsigned int GetNumberOfKernels(); std::pair< double, double > GetFaRange(){ return m_FaRange; } std::pair< double, double > GetAdcRange(){ return m_AdcRange; } unsigned int GetMaxNumKernels(){ return m_MaxNumKernels; } void Clear(); - std::vector< GradientType > m_PrototypeMaxDirection; std::vector< vnl_vector< double > > GetShCoefficients(){ return m_ShCoefficients; } std::vector< double > GetB0Signals(){ return m_B0Signal; } vnl_matrix GetSphericalCoordinates(){ return m_SphCoords; } unsigned int GetShOrder(){ return m_ShOrder; } int GetModelIndex(){ return m_ModelIndex; } double GetBaselineSignal(int index){ return m_B0Signal.at(index); } vnl_vector< double > GetCoefficients(int index){ return m_ShCoefficients.at(index); } + bool SampleKernels(DiffusionImage::Pointer diffImg, ItkUcharImageType::Pointer maskImage, TensorImageType::Pointer tensorImage=NULL, QballFilterType::CoefficientImageType::Pointer coeffImage=NULL, ItkDoubleImageType::Pointer adcImage=NULL); + protected: void Cart2Sph( GradientListType gradients ); void RandomModel(); - std::pair< double, double > m_AdcRange; - std::pair< double, double > m_FaRange; std::vector< vnl_vector< double > > m_ShCoefficients; std::vector< double > m_B0Signal; + std::vector< GradientType > m_PrototypeMaxDirection; vnl_matrix m_SphCoords; + std::pair< double, double > m_AdcRange; + std::pair< double, double > m_FaRange; unsigned int m_ShOrder; int m_ModelIndex; unsigned int m_MaxNumKernels; }; } #include "mitkRawShModel.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxSignalGenerationTest.cpp b/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxSignalGenerationTest.cpp index 1f28408a0a..6e45aa8009 100644 --- a/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxSignalGenerationTest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxSignalGenerationTest.cpp @@ -1,276 +1,277 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include +#include #include #include #include #include #include #include #include #include #include #include #include /**Documentation * Test the Fiberfox simulation functions (fiberBundle -> diffusion weighted image) */ bool CompareDwi(itk::VectorImage< short, 3 >* dwi1, itk::VectorImage< short, 3 >* dwi2) { typedef itk::VectorImage< short, 3 > DwiImageType; try{ itk::ImageRegionIterator< DwiImageType > it1(dwi1, dwi1->GetLargestPossibleRegion()); itk::ImageRegionIterator< DwiImageType > it2(dwi2, dwi2->GetLargestPossibleRegion()); while(!it1.IsAtEnd()) { if (it1.Get()!=it2.Get()) return false; ++it1; ++it2; } } catch(...) { return false; } return true; } void StartSimulation(FiberfoxParameters parameters, FiberBundleX::Pointer fiberBundle, mitk::DiffusionImage::Pointer refImage, string message) { itk::TractsToDWIImageFilter< short >::Pointer tractsToDwiFilter = itk::TractsToDWIImageFilter< short >::New(); tractsToDwiFilter->SetUseConstantRandSeed(true); tractsToDwiFilter->SetParameters(parameters); tractsToDwiFilter->SetFiberBundle(fiberBundle); tractsToDwiFilter->Update(); mitk::DiffusionImage::Pointer testImage = mitk::DiffusionImage::New(); testImage->SetVectorImage( tractsToDwiFilter->GetOutput() ); testImage->SetReferenceBValue(parameters.m_SignalGen.m_Bvalue); testImage->SetDirections(parameters.m_SignalGen.GetGradientDirections()); testImage->InitializeFromVectorImage(); if (refImage.IsNotNull()) { bool cond = CompareDwi(testImage->GetVectorImage(), refImage->GetVectorImage()); if (!cond) { MITK_INFO << "Saving test and rference image to " << mitk::IOUtil::GetTempPath(); mitk::IOUtil::SaveBaseData(testImage, mitk::IOUtil::GetTempPath()+"testImage.dwi"); mitk::IOUtil::SaveBaseData(refImage, mitk::IOUtil::GetTempPath()+"refImage.dwi"); } MITK_TEST_CONDITION_REQUIRED(cond, message); } } int mitkFiberfoxSignalGenerationTest(int argc, char* argv[]) { MITK_TEST_BEGIN("mitkFiberfoxSignalGenerationTest"); MITK_TEST_CONDITION_REQUIRED(argc>=19,"check for input data"); // input fiber bundle FiberBundleXReader::Pointer fibReader = FiberBundleXReader::New(); fibReader->SetFileName(argv[1]); fibReader->Update(); FiberBundleX::Pointer fiberBundle = dynamic_cast(fibReader->GetOutput()); // reference diffusion weighted images mitk::DiffusionImage::Pointer stickBall = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[2])->GetData()); mitk::DiffusionImage::Pointer stickAstrosticks = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[3])->GetData()); mitk::DiffusionImage::Pointer stickDot = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[4])->GetData()); mitk::DiffusionImage::Pointer tensorBall = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[5])->GetData()); mitk::DiffusionImage::Pointer stickTensorBall = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[6])->GetData()); mitk::DiffusionImage::Pointer stickTensorBallAstrosticks = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[7])->GetData()); mitk::DiffusionImage::Pointer gibbsringing = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[8])->GetData()); mitk::DiffusionImage::Pointer ghost = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[9])->GetData()); mitk::DiffusionImage::Pointer aliasing = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[10])->GetData()); mitk::DiffusionImage::Pointer eddy = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[11])->GetData()); mitk::DiffusionImage::Pointer linearmotion = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[12])->GetData()); mitk::DiffusionImage::Pointer randommotion = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[13])->GetData()); mitk::DiffusionImage::Pointer spikes = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[14])->GetData()); mitk::DiffusionImage::Pointer riciannoise = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[15])->GetData()); mitk::DiffusionImage::Pointer chisquarenoise = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[16])->GetData()); mitk::DiffusionImage::Pointer distortions = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[17])->GetData()); mitk::Image::Pointer mitkFMap = dynamic_cast(mitk::IOUtil::LoadDataNode(argv[18])->GetData()); typedef itk::Image ItkDoubleImgType; ItkDoubleImgType::Pointer fMap = ItkDoubleImgType::New(); mitk::CastToItkImage(mitkFMap, fMap); FiberfoxParameters parameters; parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; parameters.m_SignalGen.m_SignalScale = 10000; parameters.m_SignalGen.m_ImageRegion = stickBall->GetVectorImage()->GetLargestPossibleRegion(); parameters.m_SignalGen.m_ImageSpacing = stickBall->GetVectorImage()->GetSpacing(); parameters.m_SignalGen.m_ImageOrigin = stickBall->GetVectorImage()->GetOrigin(); parameters.m_SignalGen.m_ImageDirection = stickBall->GetVectorImage()->GetDirection(); parameters.m_SignalGen.m_Bvalue = stickBall->GetReferenceBValue(); parameters.m_SignalGen.SetGradienDirections(stickBall->GetDirections()); // intra and inter axonal compartments mitk::StickModel stickModel; stickModel.SetBvalue(parameters.m_SignalGen.m_Bvalue); stickModel.SetT2(110); stickModel.SetDiffusivity(0.001); stickModel.SetGradientList(parameters.m_SignalGen.GetGradientDirections()); mitk::TensorModel tensorModel; tensorModel.SetT2(110); stickModel.SetBvalue(parameters.m_SignalGen.m_Bvalue); tensorModel.SetDiffusivity1(0.001); tensorModel.SetDiffusivity2(0.00025); tensorModel.SetDiffusivity3(0.00025); tensorModel.SetGradientList(parameters.m_SignalGen.GetGradientDirections()); // extra axonal compartment models mitk::BallModel ballModel; ballModel.SetT2(80); ballModel.SetBvalue(parameters.m_SignalGen.m_Bvalue); ballModel.SetDiffusivity(0.001); ballModel.SetGradientList(parameters.m_SignalGen.GetGradientDirections()); mitk::AstroStickModel astrosticksModel; astrosticksModel.SetT2(80); astrosticksModel.SetBvalue(parameters.m_SignalGen.m_Bvalue); astrosticksModel.SetDiffusivity(0.001); astrosticksModel.SetRandomizeSticks(true); astrosticksModel.SetSeed(0); astrosticksModel.SetGradientList(parameters.m_SignalGen.GetGradientDirections()); mitk::DotModel dotModel; dotModel.SetT2(80); dotModel.SetGradientList(parameters.m_SignalGen.GetGradientDirections()); // noise models mitk::RicianNoiseModel* ricianNoiseModel = new mitk::RicianNoiseModel(); ricianNoiseModel->SetNoiseVariance(1000000); ricianNoiseModel->SetSeed(0); // Rician noise mitk::ChiSquareNoiseModel* chiSquareNoiseModel = new mitk::ChiSquareNoiseModel(); chiSquareNoiseModel->SetNoiseVariance(1000000); chiSquareNoiseModel->SetSeed(0); try{ // Stick-Ball parameters.m_FiberModelList.push_back(&stickModel); parameters.m_NonFiberModelList.push_back(&ballModel); StartSimulation(parameters, fiberBundle, stickBall, argv[2]); // Srick-Astrosticks parameters.m_NonFiberModelList.clear(); parameters.m_NonFiberModelList.push_back(&astrosticksModel); StartSimulation(parameters, fiberBundle, stickAstrosticks, argv[3]); // Stick-Dot parameters.m_NonFiberModelList.clear(); parameters.m_NonFiberModelList.push_back(&dotModel); StartSimulation(parameters, fiberBundle, stickDot, argv[4]); // Tensor-Ball parameters.m_FiberModelList.clear(); parameters.m_FiberModelList.push_back(&tensorModel); parameters.m_NonFiberModelList.clear(); parameters.m_NonFiberModelList.push_back(&ballModel); StartSimulation(parameters, fiberBundle, tensorBall, argv[5]); // Stick-Tensor-Ball parameters.m_FiberModelList.clear(); parameters.m_FiberModelList.push_back(&stickModel); parameters.m_FiberModelList.push_back(&tensorModel); parameters.m_NonFiberModelList.clear(); parameters.m_NonFiberModelList.push_back(&ballModel); StartSimulation(parameters, fiberBundle, stickTensorBall, argv[6]); // Stick-Tensor-Ball-Astrosticks parameters.m_NonFiberModelList.push_back(&astrosticksModel); StartSimulation(parameters, fiberBundle, stickTensorBallAstrosticks, argv[7]); // Gibbs ringing parameters.m_FiberModelList.clear(); parameters.m_FiberModelList.push_back(&stickModel); parameters.m_NonFiberModelList.clear(); parameters.m_NonFiberModelList.push_back(&ballModel); parameters.m_SignalGen.m_DoAddGibbsRinging = true; StartSimulation(parameters, fiberBundle, gibbsringing, argv[8]); // Ghost parameters.m_SignalGen.m_DoAddGibbsRinging = false; parameters.m_SignalGen.m_KspaceLineOffset = 0.25; StartSimulation(parameters, fiberBundle, ghost, argv[9]); // Aliasing parameters.m_SignalGen.m_KspaceLineOffset = 0; parameters.m_SignalGen.m_CroppingFactor = 0.4; parameters.m_SignalGen.m_SignalScale = 1000; StartSimulation(parameters, fiberBundle, aliasing, argv[10]); // Eddy currents parameters.m_SignalGen.m_CroppingFactor = 1; parameters.m_SignalGen.m_SignalScale = 10000; parameters.m_SignalGen.m_EddyStrength = 0.05; StartSimulation(parameters, fiberBundle, eddy, argv[11]); // Motion (linear) parameters.m_SignalGen.m_EddyStrength = 0.0; parameters.m_SignalGen.m_DoAddMotion = true; parameters.m_SignalGen.m_DoRandomizeMotion = false; parameters.m_SignalGen.m_Translation[1] = 10; parameters.m_SignalGen.m_Rotation[2] = 90; StartSimulation(parameters, fiberBundle, linearmotion, argv[12]); // Motion (random) parameters.m_SignalGen.m_DoRandomizeMotion = true; parameters.m_SignalGen.m_Translation[1] = 5; parameters.m_SignalGen.m_Rotation[2] = 45; StartSimulation(parameters, fiberBundle, randommotion, argv[13]); // Spikes parameters.m_SignalGen.m_DoAddMotion = false; parameters.m_SignalGen.m_Spikes = 5; parameters.m_SignalGen.m_SpikeAmplitude = 1; StartSimulation(parameters, fiberBundle, spikes, argv[14]); // Rician noise parameters.m_SignalGen.m_Spikes = 0; parameters.m_NoiseModel = ricianNoiseModel; StartSimulation(parameters, fiberBundle, riciannoise, argv[15]); delete parameters.m_NoiseModel; // Chi-square noise parameters.m_NoiseModel = chiSquareNoiseModel; StartSimulation(parameters, fiberBundle, chisquarenoise, argv[16]); delete parameters.m_NoiseModel; // Distortions parameters.m_NoiseModel = NULL; parameters.m_SignalGen.m_FrequencyMap = fMap; StartSimulation(parameters, fiberBundle, distortions, argv[17]); } catch (std::exception &e) { MITK_TEST_CONDITION_REQUIRED(false, e.what()); } // always end with this! MITK_TEST_END(); } diff --git a/Modules/DiffusionImaging/FiberTracking/files.cmake b/Modules/DiffusionImaging/FiberTracking/files.cmake index c0f05f474b..456eca1f99 100644 --- a/Modules/DiffusionImaging/FiberTracking/files.cmake +++ b/Modules/DiffusionImaging/FiberTracking/files.cmake @@ -1,79 +1,81 @@ set(CPP_FILES ## IO datastructures IODataStructures/FiberBundleX/mitkFiberBundleX.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.cpp IODataStructures/FiberBundleX/mitkTrackvis.cpp IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp + IODataStructures/mitkFiberfoxNoTemplateParameters.cpp # Interactions Interactions/mitkFiberBundleInteractor.cpp # Tractography Algorithms/GibbsTracking/mitkParticleGrid.cpp Algorithms/GibbsTracking/mitkMetropolisHastingsSampler.cpp Algorithms/GibbsTracking/mitkEnergyComputer.cpp Algorithms/GibbsTracking/mitkGibbsEnergyComputer.cpp Algorithms/GibbsTracking/mitkFiberBuilder.cpp Algorithms/GibbsTracking/mitkSphereInterpolator.cpp ) set(H_FILES # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.h IODataStructures/FiberBundleX/mitkFiberBundleXWriter.h IODataStructures/FiberBundleX/mitkFiberBundleXReader.h IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.h IODataStructures/FiberBundleX/mitkTrackvis.h IODataStructures/mitkFiberfoxParameters.h + IODataStructures/mitkFiberfoxNoTemplateParameters.h # Algorithms Algorithms/itkTractDensityImageFilter.h Algorithms/itkTractsToFiberEndingsImageFilter.h Algorithms/itkTractsToRgbaImageFilter.h Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h Algorithms/itkFibersFromPlanarFiguresFilter.h Algorithms/itkTractsToDWIImageFilter.h Algorithms/itkTractsToVectorImageFilter.h Algorithms/itkKspaceImageFilter.h Algorithms/itkDftImageFilter.h Algorithms/itkAddArtifactsToDwiImageFilter.h Algorithms/itkFieldmapGeneratorFilter.h Algorithms/itkEvaluateDirectionImagesFilter.h Algorithms/itkEvaluateTractogramDirectionsFilter.h # (old) Tractography Algorithms/itkGibbsTrackingFilter.h Algorithms/itkStochasticTractographyFilter.h Algorithms/itkStreamlineTrackingFilter.h Algorithms/GibbsTracking/mitkParticle.h Algorithms/GibbsTracking/mitkParticleGrid.h Algorithms/GibbsTracking/mitkMetropolisHastingsSampler.h Algorithms/GibbsTracking/mitkSimpSamp.h Algorithms/GibbsTracking/mitkEnergyComputer.h Algorithms/GibbsTracking/mitkGibbsEnergyComputer.h Algorithms/GibbsTracking/mitkSphereInterpolator.h Algorithms/GibbsTracking/mitkFiberBuilder.h # Signal Models SignalModels/mitkDiffusionSignalModel.h SignalModels/mitkTensorModel.h SignalModels/mitkBallModel.h SignalModels/mitkDotModel.h SignalModels/mitkAstroStickModel.h SignalModels/mitkStickModel.h SignalModels/mitkRawShModel.h SignalModels/mitkDiffusionNoiseModel.h SignalModels/mitkRicianNoiseModel.h SignalModels/mitkChiSquareNoiseModel.h ) set(RESOURCE_FILES # Binary directory resources FiberTrackingLUTBaryCoords.bin FiberTrackingLUTIndices.bin # Shaders Shaders/mitkShaderFiberClipping.xml ) diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.cpp index f15fc8c968..4f6c2f3b89 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.cpp @@ -1,2643 +1,2551 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ //misc #define _USE_MATH_DEFINES #include // Blueberry #include #include // Qmitk #include "QmitkFiberfoxView.h" // MITK #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "usModuleRegistry.h" #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include QmitkFiberfoxWorker::QmitkFiberfoxWorker(QmitkFiberfoxView* view) : m_View(view) { } void QmitkFiberfoxWorker::run() { try{ switch (m_FilterType) { case 0: m_View->m_TractsToDwiFilter->Update(); break; case 1: m_View->m_ArtifactsToDwiFilter->Update(); break; } } catch( ... ) { } m_View->m_Thread.quit(); } const std::string QmitkFiberfoxView::VIEW_ID = "org.mitk.views.fiberfoxview"; QmitkFiberfoxView::QmitkFiberfoxView() : QmitkAbstractView() , m_Controls( 0 ) , m_SelectedImage( NULL ) , m_Worker(this) , m_ThreadIsRunning(false) { m_Worker.moveToThread(&m_Thread); connect(&m_Thread, SIGNAL(started()), this, SLOT(BeforeThread())); connect(&m_Thread, SIGNAL(started()), &m_Worker, SLOT(run())); connect(&m_Thread, SIGNAL(finished()), this, SLOT(AfterThread())); connect(&m_Thread, SIGNAL(terminated()), this, SLOT(AfterThread())); m_SimulationTimer = new QTimer(this); } void QmitkFiberfoxView::KillThread() { MITK_INFO << "Aborting DWI simulation."; switch (m_Worker.m_FilterType) { case 0: m_TractsToDwiFilter->SetAbortGenerateData(true); break; case 1: m_ArtifactsToDwiFilter->SetAbortGenerateData(true); break; } m_Controls->m_AbortSimulationButton->setEnabled(false); m_Controls->m_AbortSimulationButton->setText("Aborting simulation ..."); } void QmitkFiberfoxView::BeforeThread() { m_SimulationTime = QTime::currentTime(); m_SimulationTimer->start(100); m_Controls->m_AbortSimulationButton->setVisible(true); m_Controls->m_GenerateImageButton->setVisible(false); m_Controls->m_SimulationStatusText->setVisible(true); m_ThreadIsRunning = true; } void QmitkFiberfoxView::AfterThread() { UpdateSimulationStatus(); m_SimulationTimer->stop(); m_Controls->m_AbortSimulationButton->setVisible(false); m_Controls->m_AbortSimulationButton->setEnabled(true); m_Controls->m_AbortSimulationButton->setText("Abort simulation"); m_Controls->m_GenerateImageButton->setVisible(true); m_ThreadIsRunning = false; QString statusText; FiberfoxParameters parameters; mitk::DiffusionImage::Pointer mitkImage = mitk::DiffusionImage::New(); switch (m_Worker.m_FilterType) { case 0: { statusText = QString(m_TractsToDwiFilter->GetStatusText().c_str()); if (m_TractsToDwiFilter->GetAbortGenerateData()) { MITK_INFO << "Simulation aborted."; return; } parameters = m_TractsToDwiFilter->GetParameters(); mitkImage->SetVectorImage( m_TractsToDwiFilter->GetOutput() ); mitkImage->SetReferenceBValue(parameters.m_SignalGen.m_Bvalue); mitkImage->SetDirections(parameters.m_SignalGen.GetGradientDirections()); mitkImage->InitializeFromVectorImage(); parameters.m_Misc.m_ResultNode->SetData( mitkImage ); parameters.m_Misc.m_ResultNode->SetName(parameters.m_Misc.m_ParentNode->GetName() +"_D"+QString::number(parameters.m_SignalGen.m_ImageRegion.GetSize(0)).toStdString() +"-"+QString::number(parameters.m_SignalGen.m_ImageRegion.GetSize(1)).toStdString() +"-"+QString::number(parameters.m_SignalGen.m_ImageRegion.GetSize(2)).toStdString() +"_S"+QString::number(parameters.m_SignalGen.m_ImageSpacing[0]).toStdString() +"-"+QString::number(parameters.m_SignalGen.m_ImageSpacing[1]).toStdString() +"-"+QString::number(parameters.m_SignalGen.m_ImageSpacing[2]).toStdString() +"_b"+QString::number(parameters.m_SignalGen.m_Bvalue).toStdString() +"_"+parameters.m_Misc.m_SignalModelString +parameters.m_Misc.m_ArtifactModelString); GetDataStorage()->Add(parameters.m_Misc.m_ResultNode, parameters.m_Misc.m_ParentNode); parameters.m_Misc.m_ResultNode->SetProperty( "levelwindow", mitk::LevelWindowProperty::New(m_TractsToDwiFilter->GetLevelWindow()) ); if (m_Controls->m_VolumeFractionsBox->isChecked()) { std::vector< itk::TractsToDWIImageFilter< short >::ItkDoubleImgType::Pointer > volumeFractions = m_TractsToDwiFilter->GetVolumeFractions(); for (unsigned int k=0; kInitializeByItk(volumeFractions.at(k).GetPointer()); image->SetVolume(volumeFractions.at(k)->GetBufferPointer()); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); node->SetName(parameters.m_Misc.m_ParentNode->GetName()+"_CompartmentVolume-"+QString::number(k).toStdString()); GetDataStorage()->Add(node, parameters.m_Misc.m_ParentNode); } } m_TractsToDwiFilter = NULL; break; } case 1: { statusText = QString(m_ArtifactsToDwiFilter->GetStatusText().c_str()); if (m_ArtifactsToDwiFilter->GetAbortGenerateData()) { MITK_INFO << "Simulation aborted."; return; } parameters = m_ArtifactsToDwiFilter->GetParameters().CopyParameters(); mitk::DiffusionImage::Pointer diffImg = dynamic_cast*>(parameters.m_Misc.m_ParentNode->GetData()); mitkImage = mitk::DiffusionImage::New(); mitkImage->SetVectorImage( m_ArtifactsToDwiFilter->GetOutput() ); mitkImage->SetReferenceBValue(diffImg->GetReferenceBValue()); mitkImage->SetDirections(diffImg->GetDirections()); mitkImage->InitializeFromVectorImage(); parameters.m_Misc.m_ResultNode->SetData( mitkImage ); parameters.m_Misc.m_ResultNode->SetName(parameters.m_Misc.m_ParentNode->GetName()+parameters.m_Misc.m_ArtifactModelString); GetDataStorage()->Add(parameters.m_Misc.m_ResultNode, parameters.m_Misc.m_ParentNode); m_ArtifactsToDwiFilter = NULL; break; } } mitk::BaseData::Pointer basedata = parameters.m_Misc.m_ResultNode->GetData(); if (basedata.IsNotNull()) { mitk::RenderingManager::GetInstance()->InitializeViews( basedata->GetTimeGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true ); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } if (!parameters.m_Misc.m_OutputPath.empty()) { try{ QString outputFileName(parameters.m_Misc.m_OutputPath.c_str()); outputFileName += parameters.m_Misc.m_ResultNode->GetName().c_str(); outputFileName.replace(QString("."), QString("_")); outputFileName += ".dwi"; QString status("Saving output image to "); status += outputFileName; m_Controls->m_SimulationStatusText->append(status); mitk::IOUtil::SaveBaseData(mitkImage, outputFileName.toStdString()); m_Controls->m_SimulationStatusText->append("File saved successfully."); } catch (itk::ExceptionObject &e) { QString status("Exception during DWI writing: "); status += e.GetDescription(); m_Controls->m_SimulationStatusText->append(status); } catch (...) { m_Controls->m_SimulationStatusText->append("Unknown exception during DWI writing!"); } } parameters.m_SignalGen.m_FrequencyMap = NULL; } void QmitkFiberfoxView::UpdateSimulationStatus() { QString statusText; switch (m_Worker.m_FilterType) { case 0: statusText = QString(m_TractsToDwiFilter->GetStatusText().c_str()); break; case 1: statusText = QString(m_ArtifactsToDwiFilter->GetStatusText().c_str()); break; } if (QString::compare(m_SimulationStatusText,statusText)!=0) { m_Controls->m_SimulationStatusText->clear(); statusText = "
"+statusText+"
"; m_Controls->m_SimulationStatusText->setText(statusText); QScrollBar *vScrollBar = m_Controls->m_SimulationStatusText->verticalScrollBar(); vScrollBar->triggerAction(QScrollBar::SliderToMaximum); } } // Destructor QmitkFiberfoxView::~QmitkFiberfoxView() { delete m_SimulationTimer; } void QmitkFiberfoxView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkFiberfoxViewControls; m_Controls->setupUi( parent ); m_Controls->m_StickWidget1->setVisible(true); m_Controls->m_StickWidget2->setVisible(false); m_Controls->m_ZeppelinWidget1->setVisible(false); m_Controls->m_ZeppelinWidget2->setVisible(false); m_Controls->m_TensorWidget1->setVisible(false); m_Controls->m_TensorWidget2->setVisible(false); m_Controls->m_BallWidget1->setVisible(true); m_Controls->m_BallWidget2->setVisible(false); m_Controls->m_AstrosticksWidget1->setVisible(false); m_Controls->m_AstrosticksWidget2->setVisible(false); m_Controls->m_DotWidget1->setVisible(false); m_Controls->m_DotWidget2->setVisible(false); m_Controls->m_PrototypeWidget1->setVisible(false); m_Controls->m_PrototypeWidget2->setVisible(false); m_Controls->m_PrototypeWidget3->setVisible(false); m_Controls->m_PrototypeWidget4->setVisible(false); m_Controls->m_PrototypeWidget3->SetMinFa(0.0); m_Controls->m_PrototypeWidget3->SetMaxFa(0.15); m_Controls->m_PrototypeWidget4->SetMinFa(0.0); m_Controls->m_PrototypeWidget4->SetMaxFa(0.15); m_Controls->m_PrototypeWidget3->SetMinAdc(0.0); m_Controls->m_PrototypeWidget3->SetMaxAdc(0.001); m_Controls->m_PrototypeWidget4->SetMinAdc(0.003); m_Controls->m_PrototypeWidget4->SetMaxAdc(0.004); m_Controls->m_Comp4FractionFrame->setVisible(false); m_Controls->m_DiffusionPropsMessage->setVisible(false); m_Controls->m_GeometryMessage->setVisible(false); m_Controls->m_AdvancedSignalOptionsFrame->setVisible(false); m_Controls->m_AdvancedFiberOptionsFrame->setVisible(false); m_Controls->m_VarianceBox->setVisible(false); m_Controls->m_NoiseFrame->setVisible(false); m_Controls->m_GhostFrame->setVisible(false); m_Controls->m_DistortionsFrame->setVisible(false); m_Controls->m_EddyFrame->setVisible(false); m_Controls->m_SpikeFrame->setVisible(false); m_Controls->m_AliasingFrame->setVisible(false); m_Controls->m_MotionArtifactFrame->setVisible(false); m_ParameterFile = QDir::currentPath()+"/param.ffp"; m_Controls->m_AbortSimulationButton->setVisible(false); m_Controls->m_SimulationStatusText->setVisible(false); m_Controls->m_FrequencyMapBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isMitkImage = mitk::TNodePredicateDataType::New(); mitk::NodePredicateDataType::Pointer isDwi = mitk::NodePredicateDataType::New("DiffusionImage"); mitk::NodePredicateDataType::Pointer isDti = mitk::NodePredicateDataType::New("TensorImage"); mitk::NodePredicateDataType::Pointer isQbi = mitk::NodePredicateDataType::New("QBallImage"); mitk::NodePredicateOr::Pointer isDiffusionImage = mitk::NodePredicateOr::New(isDwi, isDti); isDiffusionImage = mitk::NodePredicateOr::New(isDiffusionImage, isQbi); mitk::NodePredicateNot::Pointer noDiffusionImage = mitk::NodePredicateNot::New(isDiffusionImage); mitk::NodePredicateAnd::Pointer finalPredicate = mitk::NodePredicateAnd::New(isMitkImage, noDiffusionImage); m_Controls->m_FrequencyMapBox->SetPredicate(finalPredicate); m_Controls->m_Comp4VolumeFraction->SetDataStorage(this->GetDataStorage()); m_Controls->m_Comp4VolumeFraction->SetPredicate(finalPredicate); connect( m_SimulationTimer, SIGNAL(timeout()), this, SLOT(UpdateSimulationStatus()) ); connect((QObject*) m_Controls->m_AbortSimulationButton, SIGNAL(clicked()), (QObject*) this, SLOT(KillThread())); connect((QObject*) m_Controls->m_GenerateImageButton, SIGNAL(clicked()), (QObject*) this, SLOT(GenerateImage())); connect((QObject*) m_Controls->m_GenerateFibersButton, SIGNAL(clicked()), (QObject*) this, SLOT(GenerateFibers())); connect((QObject*) m_Controls->m_CircleButton, SIGNAL(clicked()), (QObject*) this, SLOT(OnDrawROI())); connect((QObject*) m_Controls->m_FlipButton, SIGNAL(clicked()), (QObject*) this, SLOT(OnFlipButton())); connect((QObject*) m_Controls->m_JoinBundlesButton, SIGNAL(clicked()), (QObject*) this, SLOT(JoinBundles())); connect((QObject*) m_Controls->m_VarianceBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnVarianceChanged(double))); connect((QObject*) m_Controls->m_DistributionBox, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(OnDistributionChanged(int))); connect((QObject*) m_Controls->m_FiberDensityBox, SIGNAL(valueChanged(int)), (QObject*) this, SLOT(OnFiberDensityChanged(int))); connect((QObject*) m_Controls->m_FiberSamplingBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnFiberSamplingChanged(double))); connect((QObject*) m_Controls->m_TensionBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnTensionChanged(double))); connect((QObject*) m_Controls->m_ContinuityBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnContinuityChanged(double))); connect((QObject*) m_Controls->m_BiasBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnBiasChanged(double))); connect((QObject*) m_Controls->m_AddNoise, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddNoise(int))); connect((QObject*) m_Controls->m_AddGhosts, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddGhosts(int))); connect((QObject*) m_Controls->m_AddDistortions, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddDistortions(int))); connect((QObject*) m_Controls->m_AddEddy, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddEddy(int))); connect((QObject*) m_Controls->m_AddSpikes, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddSpikes(int))); connect((QObject*) m_Controls->m_AddAliasing, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddAliasing(int))); connect((QObject*) m_Controls->m_AddMotion, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddMotion(int))); connect((QObject*) m_Controls->m_ConstantRadiusBox, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnConstantRadius(int))); connect((QObject*) m_Controls->m_CopyBundlesButton, SIGNAL(clicked()), (QObject*) this, SLOT(CopyBundles())); connect((QObject*) m_Controls->m_TransformBundlesButton, SIGNAL(clicked()), (QObject*) this, SLOT(ApplyTransform())); connect((QObject*) m_Controls->m_AlignOnGrid, SIGNAL(clicked()), (QObject*) this, SLOT(AlignOnGrid())); connect((QObject*) m_Controls->m_Compartment1Box, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(Comp1ModelFrameVisibility(int))); connect((QObject*) m_Controls->m_Compartment2Box, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(Comp2ModelFrameVisibility(int))); connect((QObject*) m_Controls->m_Compartment3Box, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(Comp3ModelFrameVisibility(int))); connect((QObject*) m_Controls->m_Compartment4Box, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(Comp4ModelFrameVisibility(int))); connect((QObject*) m_Controls->m_AdvancedOptionsBox, SIGNAL( stateChanged(int)), (QObject*) this, SLOT(ShowAdvancedOptions(int))); connect((QObject*) m_Controls->m_AdvancedOptionsBox_2, SIGNAL( stateChanged(int)), (QObject*) this, SLOT(ShowAdvancedOptions(int))); connect((QObject*) m_Controls->m_SaveParametersButton, SIGNAL(clicked()), (QObject*) this, SLOT(SaveParameters())); connect((QObject*) m_Controls->m_LoadParametersButton, SIGNAL(clicked()), (QObject*) this, SLOT(LoadParameters())); connect((QObject*) m_Controls->m_OutputPathButton, SIGNAL(clicked()), (QObject*) this, SLOT(SetOutputPath())); } } template< class ScalarType > FiberfoxParameters< ScalarType > QmitkFiberfoxView::UpdateImageParameters() { FiberfoxParameters< ScalarType > parameters; parameters.m_Misc.m_OutputPath = ""; parameters.m_Misc.m_CheckAdvancedFiberOptionsBox = m_Controls->m_AdvancedOptionsBox->isChecked(); parameters.m_Misc.m_CheckAdvancedSignalOptionsBox = m_Controls->m_AdvancedOptionsBox_2->isChecked(); string outputPath = m_Controls->m_SavePathEdit->text().toStdString(); if (outputPath.compare("-")!=0) { parameters.m_Misc.m_OutputPath = outputPath; parameters.m_Misc.m_OutputPath += "/"; } if (m_MaskImageNode.IsNotNull()) { mitk::Image::Pointer mitkMaskImage = dynamic_cast(m_MaskImageNode->GetData()); mitk::CastToItkImage(mitkMaskImage, parameters.m_SignalGen.m_MaskImage); itk::ImageDuplicator::Pointer duplicator = itk::ImageDuplicator::New(); duplicator->SetInputImage(parameters.m_SignalGen.m_MaskImage); duplicator->Update(); parameters.m_SignalGen.m_MaskImage = duplicator->GetOutput(); } if (m_SelectedDWI.IsNotNull()) // use parameters of selected DWI { mitk::DiffusionImage::Pointer dwi = dynamic_cast*>(m_SelectedDWI->GetData()); parameters.m_SignalGen.m_ImageRegion = dwi->GetVectorImage()->GetLargestPossibleRegion(); parameters.m_SignalGen.m_ImageSpacing = dwi->GetVectorImage()->GetSpacing(); parameters.m_SignalGen.m_ImageOrigin = dwi->GetVectorImage()->GetOrigin(); parameters.m_SignalGen.m_ImageDirection = dwi->GetVectorImage()->GetDirection(); parameters.m_SignalGen.m_Bvalue = dwi->GetReferenceBValue(); parameters.m_SignalGen.SetGradienDirections(dwi->GetDirections()); } else if (m_SelectedImage.IsNotNull()) // use geometry of selected image { mitk::Image::Pointer img = dynamic_cast(m_SelectedImage->GetData()); itk::Image< float, 3 >::Pointer itkImg = itk::Image< float, 3 >::New(); CastToItkImage< itk::Image< float, 3 > >(img, itkImg); parameters.m_SignalGen.m_ImageRegion = itkImg->GetLargestPossibleRegion(); parameters.m_SignalGen.m_ImageSpacing = itkImg->GetSpacing(); parameters.m_SignalGen.m_ImageOrigin = itkImg->GetOrigin(); parameters.m_SignalGen.m_ImageDirection = itkImg->GetDirection(); parameters.m_SignalGen.SetNumWeightedVolumes(m_Controls->m_NumGradientsBox->value()); parameters.m_SignalGen.m_Bvalue = m_Controls->m_BvalueBox->value(); } else // use GUI parameters { parameters.m_SignalGen.m_ImageRegion.SetSize(0, m_Controls->m_SizeX->value()); parameters.m_SignalGen.m_ImageRegion.SetSize(1, m_Controls->m_SizeY->value()); parameters.m_SignalGen.m_ImageRegion.SetSize(2, m_Controls->m_SizeZ->value()); parameters.m_SignalGen.m_ImageSpacing[0] = m_Controls->m_SpacingX->value(); parameters.m_SignalGen.m_ImageSpacing[1] = m_Controls->m_SpacingY->value(); parameters.m_SignalGen.m_ImageSpacing[2] = m_Controls->m_SpacingZ->value(); parameters.m_SignalGen.m_ImageOrigin[0] = parameters.m_SignalGen.m_ImageSpacing[0]/2; parameters.m_SignalGen.m_ImageOrigin[1] = parameters.m_SignalGen.m_ImageSpacing[1]/2; parameters.m_SignalGen.m_ImageOrigin[2] = parameters.m_SignalGen.m_ImageSpacing[2]/2; parameters.m_SignalGen.m_ImageDirection.SetIdentity(); parameters.m_SignalGen.SetNumWeightedVolumes(m_Controls->m_NumGradientsBox->value()); parameters.m_SignalGen.m_Bvalue = m_Controls->m_BvalueBox->value(); parameters.m_SignalGen.GenerateGradientHalfShell(); } // signal relaxation parameters.m_SignalGen.m_DoSimulateRelaxation = m_Controls->m_RelaxationBox->isChecked(); parameters.m_SignalGen.m_SimulateKspaceAcquisition = parameters.m_SignalGen.m_DoSimulateRelaxation; if (parameters.m_SignalGen.m_DoSimulateRelaxation && m_SelectedBundles.size()>0 ) parameters.m_Misc.m_ArtifactModelString += "_RELAX"; // N/2 ghosts parameters.m_Misc.m_CheckAddGhostsBox = m_Controls->m_AddGhosts->isChecked(); if (m_Controls->m_AddGhosts->isChecked()) { parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; parameters.m_Misc.m_ArtifactModelString += "_GHOST"; parameters.m_SignalGen.m_KspaceLineOffset = m_Controls->m_kOffsetBox->value(); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Ghost", DoubleProperty::New(parameters.m_SignalGen.m_KspaceLineOffset)); } else parameters.m_SignalGen.m_KspaceLineOffset = 0; // Aliasing parameters.m_Misc.m_CheckAddAliasingBox = m_Controls->m_AddAliasing->isChecked(); if (m_Controls->m_AddAliasing->isChecked()) { parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; parameters.m_Misc.m_ArtifactModelString += "_ALIASING"; parameters.m_SignalGen.m_CroppingFactor = (100-m_Controls->m_WrapBox->value())/100; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Aliasing", DoubleProperty::New(m_Controls->m_WrapBox->value())); } // Spikes parameters.m_Misc.m_CheckAddSpikesBox = m_Controls->m_AddSpikes->isChecked(); if (m_Controls->m_AddSpikes->isChecked()) { parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; parameters.m_SignalGen.m_Spikes = m_Controls->m_SpikeNumBox->value(); parameters.m_SignalGen.m_SpikeAmplitude = m_Controls->m_SpikeScaleBox->value(); parameters.m_Misc.m_ArtifactModelString += "_SPIKES"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Spikes.Number", IntProperty::New(parameters.m_SignalGen.m_Spikes)); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Spikes.Amplitude", DoubleProperty::New(parameters.m_SignalGen.m_SpikeAmplitude)); } // gibbs ringing parameters.m_SignalGen.m_DoAddGibbsRinging = m_Controls->m_AddGibbsRinging->isChecked(); if (m_Controls->m_AddGibbsRinging->isChecked()) { parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Ringing", BoolProperty::New(true)); parameters.m_Misc.m_ArtifactModelString += "_RINGING"; } // add distortions parameters.m_Misc.m_CheckAddDistortionsBox = m_Controls->m_AddDistortions->isChecked(); if (m_Controls->m_AddDistortions->isChecked() && m_Controls->m_FrequencyMapBox->GetSelectedNode().IsNotNull()) { mitk::DataNode::Pointer fMapNode = m_Controls->m_FrequencyMapBox->GetSelectedNode(); mitk::Image* img = dynamic_cast(fMapNode->GetData()); ItkDoubleImgType::Pointer itkImg = ItkDoubleImgType::New(); CastToItkImage< ItkDoubleImgType >(img, itkImg); if (m_SelectedImage.IsNull()) // use geometry of frequency map { parameters.m_SignalGen.m_ImageRegion = itkImg->GetLargestPossibleRegion(); parameters.m_SignalGen.m_ImageSpacing = itkImg->GetSpacing(); parameters.m_SignalGen.m_ImageOrigin = itkImg->GetOrigin(); parameters.m_SignalGen.m_ImageDirection = itkImg->GetDirection(); } if (parameters.m_SignalGen.m_ImageRegion.GetSize(0)==itkImg->GetLargestPossibleRegion().GetSize(0) && parameters.m_SignalGen.m_ImageRegion.GetSize(1)==itkImg->GetLargestPossibleRegion().GetSize(1) && parameters.m_SignalGen.m_ImageRegion.GetSize(2)==itkImg->GetLargestPossibleRegion().GetSize(2)) { parameters.m_SignalGen.m_SimulateKspaceAcquisition = true; itk::ImageDuplicator::Pointer duplicator = itk::ImageDuplicator::New(); duplicator->SetInputImage(itkImg); duplicator->Update(); parameters.m_SignalGen.m_FrequencyMap = duplicator->GetOutput(); parameters.m_Misc.m_ArtifactModelString += "_DISTORTED"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Distortions", BoolProperty::New(true)); } } parameters.m_SignalGen.m_EddyStrength = 0; parameters.m_Misc.m_CheckAddEddyCurrentsBox = m_Controls->m_AddEddy->isChecked(); if (m_Controls->m_AddEddy->isChecked()) { parameters.m_SignalGen.m_EddyStrength = m_Controls->m_EddyGradientStrength->value(); parameters.m_Misc.m_ArtifactModelString += "_EDDY"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Eddy-strength", DoubleProperty::New(parameters.m_SignalGen.m_EddyStrength)); } // Motion parameters.m_SignalGen.m_DoAddMotion = m_Controls->m_AddMotion->isChecked(); parameters.m_SignalGen.m_DoRandomizeMotion = m_Controls->m_RandomMotion->isChecked(); parameters.m_SignalGen.m_Translation[0] = m_Controls->m_MaxTranslationBoxX->value(); parameters.m_SignalGen.m_Translation[1] = m_Controls->m_MaxTranslationBoxY->value(); parameters.m_SignalGen.m_Translation[2] = m_Controls->m_MaxTranslationBoxZ->value(); parameters.m_SignalGen.m_Rotation[0] = m_Controls->m_MaxRotationBoxX->value(); parameters.m_SignalGen.m_Rotation[1] = m_Controls->m_MaxRotationBoxY->value(); parameters.m_SignalGen.m_Rotation[2] = m_Controls->m_MaxRotationBoxZ->value(); if ( m_Controls->m_AddMotion->isChecked() && m_SelectedBundles.size()>0 ) { parameters.m_Misc.m_ArtifactModelString += "_MOTION"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Random", BoolProperty::New(parameters.m_SignalGen.m_DoRandomizeMotion)); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Translation-x", DoubleProperty::New(parameters.m_SignalGen.m_Translation[0])); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Translation-y", DoubleProperty::New(parameters.m_SignalGen.m_Translation[1])); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Translation-z", DoubleProperty::New(parameters.m_SignalGen.m_Translation[2])); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Rotation-x", DoubleProperty::New(parameters.m_SignalGen.m_Rotation[0])); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Rotation-y", DoubleProperty::New(parameters.m_SignalGen.m_Rotation[1])); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Motion.Rotation-z", DoubleProperty::New(parameters.m_SignalGen.m_Rotation[2])); } // other imaging parameters parameters.m_SignalGen.m_tLine = m_Controls->m_LineReadoutTimeBox->value(); parameters.m_SignalGen.m_tInhom = m_Controls->m_T2starBox->value(); parameters.m_SignalGen.m_tEcho = m_Controls->m_TEbox->value(); parameters.m_SignalGen.m_DoDisablePartialVolume = m_Controls->m_EnforcePureFiberVoxelsBox->isChecked(); parameters.m_SignalGen.m_AxonRadius = m_Controls->m_FiberRadius->value(); parameters.m_SignalGen.m_SignalScale = m_Controls->m_SignalScaleBox->value(); // adjust echo time if needed if ( parameters.m_SignalGen.m_tEcho < parameters.m_SignalGen.m_ImageRegion.GetSize(1)*parameters.m_SignalGen.m_tLine ) { this->m_Controls->m_TEbox->setValue( parameters.m_SignalGen.m_ImageRegion.GetSize(1)*parameters.m_SignalGen.m_tLine ); parameters.m_SignalGen.m_tEcho = m_Controls->m_TEbox->value(); QMessageBox::information( NULL, "Warning", "Echo time is too short! Time not sufficient to read slice. Automaticall adjusted to "+QString::number(parameters.m_SignalGen.m_tEcho)+" ms"); } // Noise parameters.m_Misc.m_CheckAddNoiseBox = m_Controls->m_AddNoise->isChecked(); if (m_Controls->m_AddNoise->isChecked()) { double noiseVariance = m_Controls->m_NoiseLevel->value(); { switch (m_Controls->m_NoiseDistributionBox->currentIndex()) { case 0: { parameters.m_NoiseModel = new mitk::RicianNoiseModel(); parameters.m_Misc.m_ArtifactModelString += "_RICIAN-"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Noise-Distribution", StringProperty::New("Rician")); break; } case 1: { parameters.m_NoiseModel = new mitk::ChiSquareNoiseModel(); parameters.m_Misc.m_ArtifactModelString += "_CHISQUARED-"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Noise-Distribution", StringProperty::New("Chi-squared")); break; } default: { parameters.m_NoiseModel = new mitk::RicianNoiseModel(); parameters.m_Misc.m_ArtifactModelString += "_RICIAN-"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Noise-Distribution", StringProperty::New("Rician")); } } } parameters.m_NoiseModel->SetNoiseVariance(noiseVariance); parameters.m_Misc.m_ArtifactModelString += QString::number(noiseVariance).toStdString(); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Noise-Variance", DoubleProperty::New(noiseVariance)); } // adjusting line readout time to the adapted image size needed for the DFT unsigned int y = parameters.m_SignalGen.m_ImageRegion.GetSize(1); y += y%2; if ( y>parameters.m_SignalGen.m_ImageRegion.GetSize(1) ) parameters.m_SignalGen.m_tLine *= (double)parameters.m_SignalGen.m_ImageRegion.GetSize(1)/y; // signal models { // compartment 1 switch (m_Controls->m_Compartment1Box->currentIndex()) { case 0: { mitk::StickModel* model = new mitk::StickModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(parameters.m_SignalGen.m_Bvalue); model->SetDiffusivity(m_Controls->m_StickWidget1->GetD()); model->SetT2(m_Controls->m_StickWidget1->GetT2()); model->m_CompartmentId = 1; parameters.m_FiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Stick"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Stick") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D", DoubleProperty::New(m_Controls->m_StickWidget1->GetD()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.T2", DoubleProperty::New(model->GetT2()) ); break; } case 1: { mitk::TensorModel* model = new mitk::TensorModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(parameters.m_SignalGen.m_Bvalue); model->SetDiffusivity1(m_Controls->m_ZeppelinWidget1->GetD1()); model->SetDiffusivity2(m_Controls->m_ZeppelinWidget1->GetD2()); model->SetDiffusivity3(m_Controls->m_ZeppelinWidget1->GetD2()); model->SetT2(m_Controls->m_ZeppelinWidget1->GetT2()); model->m_CompartmentId = 1; parameters.m_FiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Zeppelin"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Zeppelin") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D1", DoubleProperty::New(m_Controls->m_ZeppelinWidget1->GetD1()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D2", DoubleProperty::New(m_Controls->m_ZeppelinWidget1->GetD2()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.T2", DoubleProperty::New(model->GetT2()) ); break; } case 2: { mitk::TensorModel* model = new mitk::TensorModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(parameters.m_SignalGen.m_Bvalue); model->SetDiffusivity1(m_Controls->m_TensorWidget1->GetD1()); model->SetDiffusivity2(m_Controls->m_TensorWidget1->GetD2()); model->SetDiffusivity3(m_Controls->m_TensorWidget1->GetD3()); model->SetT2(m_Controls->m_TensorWidget1->GetT2()); model->m_CompartmentId = 1; parameters.m_FiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Tensor"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Tensor") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D1", DoubleProperty::New(m_Controls->m_TensorWidget1->GetD1()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D2", DoubleProperty::New(m_Controls->m_TensorWidget1->GetD2()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.D3", DoubleProperty::New(m_Controls->m_TensorWidget1->GetD3()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.T2", DoubleProperty::New(model->GetT2()) ); break; } case 3: { mitk::RawShModel* model = new mitk::RawShModel(); parameters.m_SignalGen.m_SimulateKspaceAcquisition = false; model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetMaxNumKernels(m_Controls->m_PrototypeWidget1->GetNumberOfSamples()); model->SetFaRange(m_Controls->m_PrototypeWidget1->GetMinFa(), m_Controls->m_PrototypeWidget1->GetMaxFa()); model->SetAdcRange(m_Controls->m_PrototypeWidget1->GetMinAdc(), m_Controls->m_PrototypeWidget1->GetMaxAdc()); model->m_CompartmentId = 1; parameters.m_FiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Prototype"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Prototype") ); break; } } // compartment 2 switch (m_Controls->m_Compartment2Box->currentIndex()) { case 0: break; case 1: { mitk::StickModel* model = new mitk::StickModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(parameters.m_SignalGen.m_Bvalue); model->SetDiffusivity(m_Controls->m_StickWidget2->GetD()); model->SetT2(m_Controls->m_StickWidget2->GetT2()); model->m_CompartmentId = 2; parameters.m_FiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Stick"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Description", StringProperty::New("Inter-axonal compartment") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Model", StringProperty::New("Stick") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D", DoubleProperty::New(m_Controls->m_StickWidget2->GetD()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.T2", DoubleProperty::New(model->GetT2()) ); break; } case 2: { mitk::TensorModel* model = new mitk::TensorModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(parameters.m_SignalGen.m_Bvalue); model->SetDiffusivity1(m_Controls->m_ZeppelinWidget2->GetD1()); model->SetDiffusivity2(m_Controls->m_ZeppelinWidget2->GetD2()); model->SetDiffusivity3(m_Controls->m_ZeppelinWidget2->GetD2()); model->SetT2(m_Controls->m_ZeppelinWidget2->GetT2()); model->m_CompartmentId = 2; parameters.m_FiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Zeppelin"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Description", StringProperty::New("Inter-axonal compartment") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Model", StringProperty::New("Zeppelin") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D1", DoubleProperty::New(m_Controls->m_ZeppelinWidget2->GetD1()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D2", DoubleProperty::New(m_Controls->m_ZeppelinWidget2->GetD2()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.T2", DoubleProperty::New(model->GetT2()) ); break; } case 3: { mitk::TensorModel* model = new mitk::TensorModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(parameters.m_SignalGen.m_Bvalue); model->SetDiffusivity1(m_Controls->m_TensorWidget2->GetD1()); model->SetDiffusivity2(m_Controls->m_TensorWidget2->GetD2()); model->SetDiffusivity3(m_Controls->m_TensorWidget2->GetD3()); model->SetT2(m_Controls->m_TensorWidget2->GetT2()); model->m_CompartmentId = 2; parameters.m_FiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Tensor"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Description", StringProperty::New("Inter-axonal compartment") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.Model", StringProperty::New("Tensor") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D1", DoubleProperty::New(m_Controls->m_TensorWidget2->GetD1()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D2", DoubleProperty::New(m_Controls->m_TensorWidget2->GetD2()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.D3", DoubleProperty::New(m_Controls->m_TensorWidget2->GetD3()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment2.T2", DoubleProperty::New(model->GetT2()) ); break; } } // compartment 3 switch (m_Controls->m_Compartment3Box->currentIndex()) { case 0: { mitk::BallModel* model = new mitk::BallModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(parameters.m_SignalGen.m_Bvalue); model->SetDiffusivity(m_Controls->m_BallWidget1->GetD()); model->SetT2(m_Controls->m_BallWidget1->GetT2()); model->m_CompartmentId = 3; parameters.m_NonFiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Ball"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Ball") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.D", DoubleProperty::New(m_Controls->m_BallWidget1->GetD()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.T2", DoubleProperty::New(model->GetT2()) ); break; } case 1: { mitk::AstroStickModel* model = new mitk::AstroStickModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(parameters.m_SignalGen.m_Bvalue); model->SetDiffusivity(m_Controls->m_AstrosticksWidget1->GetD()); model->SetT2(m_Controls->m_AstrosticksWidget1->GetT2()); model->SetRandomizeSticks(m_Controls->m_AstrosticksWidget1->GetRandomizeSticks()); model->m_CompartmentId = 3; parameters.m_NonFiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Astrosticks"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Astrosticks") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.D", DoubleProperty::New(m_Controls->m_AstrosticksWidget1->GetD()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.T2", DoubleProperty::New(model->GetT2()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.RandomSticks", BoolProperty::New(m_Controls->m_AstrosticksWidget1->GetRandomizeSticks()) ); break; } case 2: { mitk::DotModel* model = new mitk::DotModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetT2(m_Controls->m_DotWidget1->GetT2()); model->m_CompartmentId = 3; parameters.m_NonFiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Dot"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Dot") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.T2", DoubleProperty::New(model->GetT2()) ); break; } case 3: { mitk::RawShModel* model = new mitk::RawShModel(); parameters.m_SignalGen.m_SimulateKspaceAcquisition = false; model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetMaxNumKernels(m_Controls->m_PrototypeWidget3->GetNumberOfSamples()); model->SetFaRange(m_Controls->m_PrototypeWidget3->GetMinFa(), m_Controls->m_PrototypeWidget3->GetMaxFa()); model->SetAdcRange(m_Controls->m_PrototypeWidget3->GetMinAdc(), m_Controls->m_PrototypeWidget3->GetMaxAdc()); model->m_CompartmentId = 3; parameters.m_NonFiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Prototype"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Prototype") ); break; } } // compartment 4 ItkDoubleImgType::Pointer comp4VolumeImage = NULL; ItkDoubleImgType::Pointer comp3VolumeImage = NULL; if (m_Controls->m_Compartment4Box->currentIndex()>0) { mitk::DataNode::Pointer volumeNode = m_Controls->m_Comp4VolumeFraction->GetSelectedNode(); if (volumeNode.IsNull()) { QMessageBox::information( NULL, "Information", "No volume fraction image selected! Second extra-axonal compartment has been disabled for this simultation."); MITK_WARN << "No volume fraction image selected! Second extra-axonal compartment has been disabled."; } else { MITK_INFO << "Rescaling volume fraction image..."; comp4VolumeImage = ItkDoubleImgType::New(); mitk::Image* img = dynamic_cast(volumeNode->GetData()); CastToItkImage< ItkDoubleImgType >(img, comp4VolumeImage); double max = itk::NumericTraits::min(); double min = itk::NumericTraits::max(); itk::ImageRegionIterator< ItkDoubleImgType > it(comp4VolumeImage, comp4VolumeImage->GetLargestPossibleRegion()); while(!it.IsAtEnd()) { if (parameters.m_SignalGen.m_MaskImage.IsNotNull() && parameters.m_SignalGen.m_MaskImage->GetPixel(it.GetIndex())<=0) { it.Set(0.0); ++it; continue; } if (it.Get()>max) max = it.Get(); if (it.Get()::Pointer scaler = itk::ShiftScaleImageFilter< ItkDoubleImgType, ItkDoubleImgType >::New(); scaler->SetInput(comp4VolumeImage); scaler->SetShift(-min); scaler->SetScale(1.0/(max-min)); scaler->Update(); comp4VolumeImage = scaler->GetOutput(); // itk::ImageFileWriter< ItkDoubleImgType >::Pointer wr = itk::ImageFileWriter< ItkDoubleImgType >::New(); // wr->SetInput(comp4VolumeImage); // wr->SetFileName("/local/comp4.nrrd"); // wr->Update(); // if (max>1 || min<0) // are volume fractions between 0 and 1? // { // itk::RescaleIntensityImageFilter::Pointer rescaler = itk::RescaleIntensityImageFilter::New(); // rescaler->SetInput(0, comp4VolumeImage); // rescaler->SetOutputMaximum(1); // rescaler->SetOutputMinimum(0); // rescaler->Update(); // comp4VolumeImage = rescaler->GetOutput(); // } itk::InvertIntensityImageFilter< ItkDoubleImgType, ItkDoubleImgType >::Pointer inverter = itk::InvertIntensityImageFilter< ItkDoubleImgType, ItkDoubleImgType >::New(); inverter->SetMaximum(1.0); inverter->SetInput(comp4VolumeImage); inverter->Update(); comp3VolumeImage = inverter->GetOutput(); } } if (comp4VolumeImage.IsNotNull()) { switch (m_Controls->m_Compartment4Box->currentIndex()) { case 0: break; case 1: { mitk::BallModel* model = new mitk::BallModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(parameters.m_SignalGen.m_Bvalue); model->SetDiffusivity(m_Controls->m_BallWidget2->GetD()); model->SetT2(m_Controls->m_BallWidget2->GetT2()); model->SetVolumeFractionImage(comp4VolumeImage); model->m_CompartmentId = 4; parameters.m_NonFiberModelList.back()->SetVolumeFractionImage(comp3VolumeImage); parameters.m_NonFiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Ball"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Ball") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.D", DoubleProperty::New(m_Controls->m_BallWidget2->GetD()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.T2", DoubleProperty::New(model->GetT2()) ); break; } case 2: { mitk::AstroStickModel* model = new mitk::AstroStickModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetBvalue(parameters.m_SignalGen.m_Bvalue); model->SetDiffusivity(m_Controls->m_AstrosticksWidget2->GetD()); model->SetT2(m_Controls->m_AstrosticksWidget2->GetT2()); model->SetRandomizeSticks(m_Controls->m_AstrosticksWidget2->GetRandomizeSticks()); parameters.m_NonFiberModelList.back()->SetVolumeFractionImage(comp3VolumeImage); model->SetVolumeFractionImage(comp4VolumeImage); model->m_CompartmentId = 4; parameters.m_NonFiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Astrosticks"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Astrosticks") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.D", DoubleProperty::New(m_Controls->m_AstrosticksWidget2->GetD()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.T2", DoubleProperty::New(model->GetT2()) ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.RandomSticks", BoolProperty::New(m_Controls->m_AstrosticksWidget2->GetRandomizeSticks()) ); break; } case 3: { mitk::DotModel* model = new mitk::DotModel(); model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetT2(m_Controls->m_DotWidget2->GetT2()); model->SetVolumeFractionImage(comp4VolumeImage); model->m_CompartmentId = 4; parameters.m_NonFiberModelList.back()->SetVolumeFractionImage(comp3VolumeImage); parameters.m_NonFiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Dot"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Dot") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.T2", DoubleProperty::New(model->GetT2()) ); break; } case 4: { mitk::RawShModel* model = new mitk::RawShModel(); parameters.m_SignalGen.m_SimulateKspaceAcquisition = false; model->SetGradientList(parameters.m_SignalGen.GetGradientDirections()); model->SetMaxNumKernels(m_Controls->m_PrototypeWidget4->GetNumberOfSamples()); model->SetFaRange(m_Controls->m_PrototypeWidget4->GetMinFa(), m_Controls->m_PrototypeWidget4->GetMaxFa()); model->SetAdcRange(m_Controls->m_PrototypeWidget4->GetMinAdc(), m_Controls->m_PrototypeWidget4->GetMaxAdc()); model->SetVolumeFractionImage(comp4VolumeImage); model->m_CompartmentId = 4; parameters.m_NonFiberModelList.back()->SetVolumeFractionImage(comp3VolumeImage); parameters.m_NonFiberModelList.push_back(model); parameters.m_Misc.m_SignalModelString += "Prototype"; parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Prototype") ); break; } } } } parameters.m_SignalGen.m_FiberSeparationThreshold = m_Controls->m_SeparationAngleBox->value(); switch (m_Controls->m_DiffusionDirectionBox->currentIndex()) { case 0: parameters.m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::FIBER_TANGENT_DIRECTIONS; break; case 1: parameters.m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::MAIN_FIBER_DIRECTIONS; break; case 2: parameters.m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::RANDOM_DIRECTIONS; parameters.m_SignalGen.m_DoAddMotion = false; parameters.m_SignalGen.m_DoAddGibbsRinging = false; parameters.m_SignalGen.m_KspaceLineOffset = 0.0; parameters.m_SignalGen.m_FrequencyMap = NULL; parameters.m_SignalGen.m_CroppingFactor = 1.0; parameters.m_SignalGen.m_EddyStrength = 0; break; default: parameters.m_SignalGen.m_DiffusionDirectionMode = SignalGenerationParameters::FIBER_TANGENT_DIRECTIONS; } parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.SignalScale", IntProperty::New(parameters.m_SignalGen.m_SignalScale)); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.FiberRadius", IntProperty::New(parameters.m_SignalGen.m_AxonRadius)); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Tinhom", DoubleProperty::New(parameters.m_SignalGen.m_tInhom)); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Tline", DoubleProperty::New(parameters.m_SignalGen.m_tLine)); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.TE", DoubleProperty::New(parameters.m_SignalGen.m_tEcho)); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.b-value", DoubleProperty::New(parameters.m_SignalGen.m_Bvalue)); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.NoPartialVolume", BoolProperty::New(parameters.m_SignalGen.m_DoDisablePartialVolume)); parameters.m_Misc.m_ResultNode->AddProperty("Fiberfox.Relaxation", BoolProperty::New(parameters.m_SignalGen.m_DoSimulateRelaxation)); parameters.m_Misc.m_ResultNode->AddProperty("binary", BoolProperty::New(false)); parameters.m_Misc.m_CheckRealTimeFibersBox = m_Controls->m_RealTimeFibers->isChecked(); parameters.m_Misc.m_CheckAdvancedFiberOptionsBox = m_Controls->m_AdvancedOptionsBox->isChecked(); parameters.m_Misc.m_CheckIncludeFiducialsBox = m_Controls->m_IncludeFiducials->isChecked(); parameters.m_Misc.m_CheckConstantRadiusBox = m_Controls->m_ConstantRadiusBox->isChecked(); switch(m_Controls->m_DistributionBox->currentIndex()) { case 0: parameters.m_FiberGen.m_Distribution = FiberGenerationParameters::DISTRIBUTE_UNIFORM; break; case 1: parameters.m_FiberGen.m_Distribution = FiberGenerationParameters::DISTRIBUTE_GAUSSIAN; break; default: parameters.m_FiberGen.m_Distribution = FiberGenerationParameters::DISTRIBUTE_UNIFORM; } parameters.m_FiberGen.m_Variance = m_Controls->m_VarianceBox->value(); parameters.m_FiberGen.m_Density = m_Controls->m_FiberDensityBox->value(); parameters.m_FiberGen.m_Sampling = m_Controls->m_FiberSamplingBox->value(); parameters.m_FiberGen.m_Tension = m_Controls->m_TensionBox->value(); parameters.m_FiberGen.m_Continuity = m_Controls->m_ContinuityBox->value(); parameters.m_FiberGen.m_Bias = m_Controls->m_BiasBox->value(); parameters.m_FiberGen.m_Rotation[0] = m_Controls->m_XrotBox->value(); parameters.m_FiberGen.m_Rotation[1] = m_Controls->m_YrotBox->value(); parameters.m_FiberGen.m_Rotation[2] = m_Controls->m_ZrotBox->value(); parameters.m_FiberGen.m_Translation[0] = m_Controls->m_XtransBox->value(); parameters.m_FiberGen.m_Translation[1] = m_Controls->m_YtransBox->value(); parameters.m_FiberGen.m_Translation[2] = m_Controls->m_ZtransBox->value(); parameters.m_FiberGen.m_Scale[0] = m_Controls->m_XscaleBox->value(); parameters.m_FiberGen.m_Scale[1] = m_Controls->m_YscaleBox->value(); parameters.m_FiberGen.m_Scale[2] = m_Controls->m_ZscaleBox->value(); return parameters; } void QmitkFiberfoxView::SaveParameters() { FiberfoxParameters<> ffParamaters = UpdateImageParameters(); QString filename = QFileDialog::getSaveFileName( 0, tr("Save Parameters"), m_ParameterFile, tr("Fiberfox Parameters (*.ffp)") ); ffParamaters.SaveParameters(filename.toStdString()); m_ParameterFile = filename; } void QmitkFiberfoxView::LoadParameters() { QString filename = QFileDialog::getOpenFileName(0, tr("Load Parameters"), QString(itksys::SystemTools::GetFilenamePath(m_ParameterFile.toStdString()).c_str()), tr("Fiberfox Parameters (*.ffp)") ); if(filename.isEmpty() || filename.isNull()) return; m_ParameterFile = filename; FiberfoxParameters<> parameters; parameters.LoadParameters(filename.toStdString()); m_Controls->m_RealTimeFibers->setChecked(parameters.m_Misc.m_CheckRealTimeFibersBox); m_Controls->m_AdvancedOptionsBox->setChecked(parameters.m_Misc.m_CheckAdvancedFiberOptionsBox); m_Controls->m_IncludeFiducials->setChecked(parameters.m_Misc.m_CheckIncludeFiducialsBox); m_Controls->m_ConstantRadiusBox->setChecked(parameters.m_Misc.m_CheckConstantRadiusBox); m_Controls->m_DistributionBox->setCurrentIndex(parameters.m_FiberGen.m_Distribution); m_Controls->m_VarianceBox->setValue(parameters.m_FiberGen.m_Variance); m_Controls->m_FiberDensityBox->setValue(parameters.m_FiberGen.m_Density); m_Controls->m_FiberSamplingBox->setValue(parameters.m_FiberGen.m_Sampling); m_Controls->m_TensionBox->setValue(parameters.m_FiberGen.m_Tension); m_Controls->m_ContinuityBox->setValue(parameters.m_FiberGen.m_Continuity); m_Controls->m_BiasBox->setValue(parameters.m_FiberGen.m_Bias); m_Controls->m_XrotBox->setValue(parameters.m_FiberGen.m_Rotation[0]); m_Controls->m_YrotBox->setValue(parameters.m_FiberGen.m_Rotation[1]); m_Controls->m_ZrotBox->setValue(parameters.m_FiberGen.m_Rotation[2]); m_Controls->m_XtransBox->setValue(parameters.m_FiberGen.m_Translation[0]); m_Controls->m_YtransBox->setValue(parameters.m_FiberGen.m_Translation[1]); m_Controls->m_ZtransBox->setValue(parameters.m_FiberGen.m_Translation[2]); m_Controls->m_XscaleBox->setValue(parameters.m_FiberGen.m_Scale[0]); m_Controls->m_YscaleBox->setValue(parameters.m_FiberGen.m_Scale[1]); m_Controls->m_ZscaleBox->setValue(parameters.m_FiberGen.m_Scale[2]); // image generation parameters m_Controls->m_SizeX->setValue(parameters.m_SignalGen.m_ImageRegion.GetSize(0)); m_Controls->m_SizeY->setValue(parameters.m_SignalGen.m_ImageRegion.GetSize(1)); m_Controls->m_SizeZ->setValue(parameters.m_SignalGen.m_ImageRegion.GetSize(2)); m_Controls->m_SpacingX->setValue(parameters.m_SignalGen.m_ImageSpacing[0]); m_Controls->m_SpacingY->setValue(parameters.m_SignalGen.m_ImageSpacing[1]); m_Controls->m_SpacingZ->setValue(parameters.m_SignalGen.m_ImageSpacing[2]); m_Controls->m_NumGradientsBox->setValue(parameters.m_SignalGen.GetNumWeightedVolumes()); m_Controls->m_BvalueBox->setValue(parameters.m_SignalGen.m_Bvalue); m_Controls->m_SignalScaleBox->setValue(parameters.m_SignalGen.m_SignalScale); m_Controls->m_TEbox->setValue(parameters.m_SignalGen.m_tEcho); m_Controls->m_LineReadoutTimeBox->setValue(parameters.m_SignalGen.m_tLine); m_Controls->m_T2starBox->setValue(parameters.m_SignalGen.m_tInhom); m_Controls->m_FiberRadius->setValue(parameters.m_SignalGen.m_AxonRadius); m_Controls->m_RelaxationBox->setChecked(parameters.m_SignalGen.m_DoSimulateRelaxation); m_Controls->m_EnforcePureFiberVoxelsBox->setChecked(parameters.m_SignalGen.m_DoDisablePartialVolume); if (parameters.m_NoiseModel!=NULL) { m_Controls->m_AddNoise->setChecked(parameters.m_Misc.m_CheckAddNoiseBox); if (dynamic_cast*>(parameters.m_NoiseModel)) m_Controls->m_NoiseDistributionBox->setCurrentIndex(0); else if (dynamic_cast*>(parameters.m_NoiseModel)) m_Controls->m_NoiseDistributionBox->setCurrentIndex(1); m_Controls->m_NoiseLevel->setValue(parameters.m_NoiseModel->GetNoiseVariance()); } else m_Controls->m_AddNoise->setChecked(false); m_Controls->m_VolumeFractionsBox->setChecked(parameters.m_Misc.m_CheckOutputVolumeFractionsBox); m_Controls->m_AdvancedOptionsBox_2->setChecked(parameters.m_Misc.m_CheckAdvancedSignalOptionsBox); m_Controls->m_AddGhosts->setChecked(parameters.m_Misc.m_CheckAddGhostsBox); m_Controls->m_AddAliasing->setChecked(parameters.m_Misc.m_CheckAddAliasingBox); m_Controls->m_AddDistortions->setChecked(parameters.m_Misc.m_CheckAddDistortionsBox); m_Controls->m_AddSpikes->setChecked(parameters.m_Misc.m_CheckAddSpikesBox); m_Controls->m_AddEddy->setChecked(parameters.m_Misc.m_CheckAddEddyCurrentsBox); m_Controls->m_kOffsetBox->setValue(parameters.m_SignalGen.m_KspaceLineOffset); m_Controls->m_WrapBox->setValue(100*(1-parameters.m_SignalGen.m_CroppingFactor)); m_Controls->m_SpikeNumBox->setValue(parameters.m_SignalGen.m_Spikes); m_Controls->m_SpikeScaleBox->setValue(parameters.m_SignalGen.m_SpikeAmplitude); m_Controls->m_EddyGradientStrength->setValue(parameters.m_SignalGen.m_EddyStrength); m_Controls->m_AddGibbsRinging->setChecked(parameters.m_SignalGen.m_DoAddGibbsRinging); m_Controls->m_AddMotion->setChecked(parameters.m_SignalGen.m_DoAddMotion); m_Controls->m_RandomMotion->setChecked(parameters.m_SignalGen.m_DoRandomizeMotion); m_Controls->m_MaxTranslationBoxX->setValue(parameters.m_SignalGen.m_Translation[0]); m_Controls->m_MaxTranslationBoxY->setValue(parameters.m_SignalGen.m_Translation[1]); m_Controls->m_MaxTranslationBoxZ->setValue(parameters.m_SignalGen.m_Translation[2]); m_Controls->m_MaxRotationBoxX->setValue(parameters.m_SignalGen.m_Rotation[0]); m_Controls->m_MaxRotationBoxY->setValue(parameters.m_SignalGen.m_Rotation[1]); m_Controls->m_MaxRotationBoxZ->setValue(parameters.m_SignalGen.m_Rotation[2]); m_Controls->m_DiffusionDirectionBox->setCurrentIndex(parameters.m_SignalGen.m_DiffusionDirectionMode); m_Controls->m_SeparationAngleBox->setValue(parameters.m_SignalGen.m_FiberSeparationThreshold); m_Controls->m_Compartment1Box->setCurrentIndex(0); m_Controls->m_Compartment2Box->setCurrentIndex(0); m_Controls->m_Compartment3Box->setCurrentIndex(0); m_Controls->m_Compartment4Box->setCurrentIndex(0); for (unsigned int i=0; i* signalModel = NULL; if (im_CompartmentId) { case 1: { if (dynamic_cast*>(signalModel)) { mitk::StickModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_StickWidget1->SetT2(model->GetT2()); m_Controls->m_StickWidget1->SetD(model->GetDiffusivity()); m_Controls->m_Compartment1Box->setCurrentIndex(0); break; } else if (dynamic_cast*>(signalModel)) { mitk::TensorModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_TensorWidget1->SetT2(model->GetT2()); m_Controls->m_TensorWidget1->SetD1(model->GetDiffusivity1()); m_Controls->m_TensorWidget1->SetD2(model->GetDiffusivity1()); m_Controls->m_TensorWidget1->SetD3(model->GetDiffusivity1()); m_Controls->m_Compartment1Box->setCurrentIndex(2); break; } else if (dynamic_cast*>(signalModel)) { mitk::RawShModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_PrototypeWidget1->SetNumberOfSamples(model->GetMaxNumKernels()); m_Controls->m_PrototypeWidget1->SetMinFa(model->GetFaRange().first); m_Controls->m_PrototypeWidget1->SetMaxFa(model->GetFaRange().second); m_Controls->m_PrototypeWidget1->SetMinAdc(model->GetAdcRange().first); m_Controls->m_PrototypeWidget1->SetMaxAdc(model->GetAdcRange().second); m_Controls->m_Compartment1Box->setCurrentIndex(3); break; } break; } case 2: { if (dynamic_cast*>(signalModel)) { mitk::StickModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_StickWidget2->SetT2(model->GetT2()); m_Controls->m_StickWidget2->SetD(model->GetDiffusivity()); m_Controls->m_Compartment2Box->setCurrentIndex(1); break; } else if (dynamic_cast*>(signalModel)) { mitk::TensorModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_TensorWidget2->SetT2(model->GetT2()); m_Controls->m_TensorWidget2->SetD1(model->GetDiffusivity1()); m_Controls->m_TensorWidget2->SetD2(model->GetDiffusivity1()); m_Controls->m_TensorWidget2->SetD3(model->GetDiffusivity1()); m_Controls->m_Compartment2Box->setCurrentIndex(3); break; } break; } case 3: { if (dynamic_cast*>(signalModel)) { mitk::BallModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_BallWidget1->SetT2(model->GetT2()); m_Controls->m_BallWidget1->SetD(model->GetDiffusivity()); m_Controls->m_Compartment3Box->setCurrentIndex(0); break; } else if (dynamic_cast*>(signalModel)) { mitk::AstroStickModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_AstrosticksWidget1->SetT2(model->GetT2()); m_Controls->m_AstrosticksWidget1->SetD(model->GetDiffusivity()); m_Controls->m_AstrosticksWidget1->SetRandomizeSticks(model->GetRandomizeSticks()); m_Controls->m_Compartment3Box->setCurrentIndex(1); break; } else if (dynamic_cast*>(signalModel)) { mitk::DotModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_DotWidget1->SetT2(model->GetT2()); m_Controls->m_Compartment3Box->setCurrentIndex(2); break; } else if (dynamic_cast*>(signalModel)) { mitk::RawShModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_PrototypeWidget3->SetNumberOfSamples(model->GetMaxNumKernels()); m_Controls->m_PrototypeWidget3->SetMinFa(model->GetFaRange().first); m_Controls->m_PrototypeWidget3->SetMaxFa(model->GetFaRange().second); m_Controls->m_PrototypeWidget3->SetMinAdc(model->GetAdcRange().first); m_Controls->m_PrototypeWidget3->SetMaxAdc(model->GetAdcRange().second); m_Controls->m_Compartment3Box->setCurrentIndex(3); break; } break; } case 4: { if (dynamic_cast*>(signalModel)) { mitk::BallModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_BallWidget2->SetT2(model->GetT2()); m_Controls->m_BallWidget2->SetD(model->GetDiffusivity()); m_Controls->m_Compartment4Box->setCurrentIndex(1); break; } else if (dynamic_cast*>(signalModel)) { mitk::AstroStickModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_AstrosticksWidget2->SetT2(model->GetT2()); m_Controls->m_AstrosticksWidget2->SetD(model->GetDiffusivity()); m_Controls->m_AstrosticksWidget2->SetRandomizeSticks(model->GetRandomizeSticks()); m_Controls->m_Compartment4Box->setCurrentIndex(2); break; } else if (dynamic_cast*>(signalModel)) { mitk::DotModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_DotWidget2->SetT2(model->GetT2()); m_Controls->m_Compartment4Box->setCurrentIndex(3); break; } else if (dynamic_cast*>(signalModel)) { mitk::RawShModel<>* model = dynamic_cast*>(signalModel); m_Controls->m_PrototypeWidget4->SetNumberOfSamples(model->GetMaxNumKernels()); m_Controls->m_PrototypeWidget4->SetMinFa(model->GetFaRange().first); m_Controls->m_PrototypeWidget4->SetMaxFa(model->GetFaRange().second); m_Controls->m_PrototypeWidget4->SetMinAdc(model->GetAdcRange().first); m_Controls->m_PrototypeWidget4->SetMaxAdc(model->GetAdcRange().second); m_Controls->m_Compartment4Box->setCurrentIndex(4); break; } break; } } } } void QmitkFiberfoxView::ShowAdvancedOptions(int state) { if (state) { m_Controls->m_AdvancedFiberOptionsFrame->setVisible(true); m_Controls->m_AdvancedSignalOptionsFrame->setVisible(true); m_Controls->m_AdvancedOptionsBox->setChecked(true); m_Controls->m_AdvancedOptionsBox_2->setChecked(true); } else { m_Controls->m_AdvancedFiberOptionsFrame->setVisible(false); m_Controls->m_AdvancedSignalOptionsFrame->setVisible(false); m_Controls->m_AdvancedOptionsBox->setChecked(false); m_Controls->m_AdvancedOptionsBox_2->setChecked(false); } } void QmitkFiberfoxView::Comp1ModelFrameVisibility(int index) { m_Controls->m_StickWidget1->setVisible(false); m_Controls->m_ZeppelinWidget1->setVisible(false); m_Controls->m_TensorWidget1->setVisible(false); m_Controls->m_PrototypeWidget1->setVisible(false); switch (index) { case 0: m_Controls->m_StickWidget1->setVisible(true); break; case 1: m_Controls->m_ZeppelinWidget1->setVisible(true); break; case 2: m_Controls->m_TensorWidget1->setVisible(true); break; case 3: m_Controls->m_PrototypeWidget1->setVisible(true); break; } } void QmitkFiberfoxView::Comp2ModelFrameVisibility(int index) { m_Controls->m_StickWidget2->setVisible(false); m_Controls->m_ZeppelinWidget2->setVisible(false); m_Controls->m_TensorWidget2->setVisible(false); switch (index) { case 0: break; case 1: m_Controls->m_StickWidget2->setVisible(true); break; case 2: m_Controls->m_ZeppelinWidget2->setVisible(true); break; case 3: m_Controls->m_TensorWidget2->setVisible(true); break; } } void QmitkFiberfoxView::Comp3ModelFrameVisibility(int index) { m_Controls->m_BallWidget1->setVisible(false); m_Controls->m_AstrosticksWidget1->setVisible(false); m_Controls->m_DotWidget1->setVisible(false); m_Controls->m_PrototypeWidget3->setVisible(false); switch (index) { case 0: m_Controls->m_BallWidget1->setVisible(true); break; case 1: m_Controls->m_AstrosticksWidget1->setVisible(true); break; case 2: m_Controls->m_DotWidget1->setVisible(true); break; case 3: m_Controls->m_PrototypeWidget3->setVisible(true); break; } } void QmitkFiberfoxView::Comp4ModelFrameVisibility(int index) { m_Controls->m_BallWidget2->setVisible(false); m_Controls->m_AstrosticksWidget2->setVisible(false); m_Controls->m_DotWidget2->setVisible(false); m_Controls->m_PrototypeWidget4->setVisible(false); m_Controls->m_Comp4FractionFrame->setVisible(false); switch (index) { case 0: break; case 1: m_Controls->m_BallWidget2->setVisible(true); m_Controls->m_Comp4FractionFrame->setVisible(true); break; case 2: m_Controls->m_AstrosticksWidget2->setVisible(true); m_Controls->m_Comp4FractionFrame->setVisible(true); break; case 3: m_Controls->m_DotWidget2->setVisible(true); m_Controls->m_Comp4FractionFrame->setVisible(true); break; case 4: m_Controls->m_PrototypeWidget4->setVisible(true); m_Controls->m_Comp4FractionFrame->setVisible(true); break; } } void QmitkFiberfoxView::OnConstantRadius(int value) { if (value>0 && m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnAddMotion(int value) { if (value>0) m_Controls->m_MotionArtifactFrame->setVisible(true); else m_Controls->m_MotionArtifactFrame->setVisible(false); } void QmitkFiberfoxView::OnAddAliasing(int value) { if (value>0) m_Controls->m_AliasingFrame->setVisible(true); else m_Controls->m_AliasingFrame->setVisible(false); } void QmitkFiberfoxView::OnAddSpikes(int value) { if (value>0) m_Controls->m_SpikeFrame->setVisible(true); else m_Controls->m_SpikeFrame->setVisible(false); } void QmitkFiberfoxView::OnAddEddy(int value) { if (value>0) m_Controls->m_EddyFrame->setVisible(true); else m_Controls->m_EddyFrame->setVisible(false); } void QmitkFiberfoxView::OnAddDistortions(int value) { if (value>0) m_Controls->m_DistortionsFrame->setVisible(true); else m_Controls->m_DistortionsFrame->setVisible(false); } void QmitkFiberfoxView::OnAddGhosts(int value) { if (value>0) m_Controls->m_GhostFrame->setVisible(true); else m_Controls->m_GhostFrame->setVisible(false); } void QmitkFiberfoxView::OnAddNoise(int value) { if (value>0) m_Controls->m_NoiseFrame->setVisible(true); else m_Controls->m_NoiseFrame->setVisible(false); } void QmitkFiberfoxView::OnDistributionChanged(int value) { if (value==1) m_Controls->m_VarianceBox->setVisible(true); else m_Controls->m_VarianceBox->setVisible(false); if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnVarianceChanged(double) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnFiberDensityChanged(int) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnFiberSamplingChanged(double) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnTensionChanged(double) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnContinuityChanged(double) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnBiasChanged(double) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::AlignOnGrid() { for (unsigned int i=0; i(m_SelectedFiducials.at(i)->GetData()); mitk::Point3D wc0 = pe->GetWorldControlPoint(0); mitk::DataStorage::SetOfObjects::ConstPointer parentFibs = GetDataStorage()->GetSources(m_SelectedFiducials.at(i)); for( mitk::DataStorage::SetOfObjects::const_iterator it = parentFibs->begin(); it != parentFibs->end(); ++it ) { mitk::DataNode::Pointer pFibNode = *it; if ( pFibNode.IsNotNull() && dynamic_cast(pFibNode->GetData()) ) { mitk::DataStorage::SetOfObjects::ConstPointer parentImgs = GetDataStorage()->GetSources(pFibNode); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = parentImgs->begin(); it2 != parentImgs->end(); ++it2 ) { mitk::DataNode::Pointer pImgNode = *it2; if ( pImgNode.IsNotNull() && dynamic_cast(pImgNode->GetData()) ) { mitk::Image::Pointer img = dynamic_cast(pImgNode->GetData()); mitk::BaseGeometry::Pointer geom = img->GetGeometry(); itk::Index<3> idx; geom->WorldToIndex(wc0, idx); mitk::Point3D cIdx; cIdx[0]=idx[0]; cIdx[1]=idx[1]; cIdx[2]=idx[2]; mitk::Point3D world; geom->IndexToWorld(cIdx,world); mitk::Vector3D trans = world - wc0; pe->GetGeometry()->Translate(trans); break; } } break; } } } for(unsigned int i=0; iGetSources(fibNode); for( mitk::DataStorage::SetOfObjects::const_iterator it = sources->begin(); it != sources->end(); ++it ) { mitk::DataNode::Pointer imgNode = *it; if ( imgNode.IsNotNull() && dynamic_cast(imgNode->GetData()) ) { mitk::DataStorage::SetOfObjects::ConstPointer derivations = GetDataStorage()->GetDerivations(fibNode); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = derivations->begin(); it2 != derivations->end(); ++it2 ) { mitk::DataNode::Pointer fiducialNode = *it2; if ( fiducialNode.IsNotNull() && dynamic_cast(fiducialNode->GetData()) ) { mitk::PlanarEllipse::Pointer pe = dynamic_cast(fiducialNode->GetData()); mitk::Point3D wc0 = pe->GetWorldControlPoint(0); mitk::Image::Pointer img = dynamic_cast(imgNode->GetData()); mitk::BaseGeometry::Pointer geom = img->GetGeometry(); itk::Index<3> idx; geom->WorldToIndex(wc0, idx); mitk::Point3D cIdx; cIdx[0]=idx[0]; cIdx[1]=idx[1]; cIdx[2]=idx[2]; mitk::Point3D world; geom->IndexToWorld(cIdx,world); mitk::Vector3D trans = world - wc0; pe->GetGeometry()->Translate(trans); } } break; } } } for(unsigned int i=0; i(m_SelectedImages.at(i)->GetData()); mitk::DataStorage::SetOfObjects::ConstPointer derivations = GetDataStorage()->GetDerivations(m_SelectedImages.at(i)); for( mitk::DataStorage::SetOfObjects::const_iterator it = derivations->begin(); it != derivations->end(); ++it ) { mitk::DataNode::Pointer fibNode = *it; if ( fibNode.IsNotNull() && dynamic_cast(fibNode->GetData()) ) { mitk::DataStorage::SetOfObjects::ConstPointer derivations2 = GetDataStorage()->GetDerivations(fibNode); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = derivations2->begin(); it2 != derivations2->end(); ++it2 ) { mitk::DataNode::Pointer fiducialNode = *it2; if ( fiducialNode.IsNotNull() && dynamic_cast(fiducialNode->GetData()) ) { mitk::PlanarEllipse::Pointer pe = dynamic_cast(fiducialNode->GetData()); mitk::Point3D wc0 = pe->GetWorldControlPoint(0); mitk::BaseGeometry::Pointer geom = img->GetGeometry(); itk::Index<3> idx; geom->WorldToIndex(wc0, idx); mitk::Point3D cIdx; cIdx[0]=idx[0]; cIdx[1]=idx[1]; cIdx[2]=idx[2]; mitk::Point3D world; geom->IndexToWorld(cIdx,world); mitk::Vector3D trans = world - wc0; pe->GetGeometry()->Translate(trans); } } } } } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnFlipButton() { if (m_SelectedFiducial.IsNull()) return; std::map::iterator it = m_DataNodeToPlanarFigureData.find(m_SelectedFiducial.GetPointer()); if( it != m_DataNodeToPlanarFigureData.end() ) { QmitkPlanarFigureData& data = it->second; data.m_Flipped += 1; data.m_Flipped %= 2; } if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } QmitkFiberfoxView::GradientListType QmitkFiberfoxView::GenerateHalfShell(int NPoints) { NPoints *= 2; GradientListType pointshell; int numB0 = NPoints/20; if (numB0==0) numB0=1; GradientType g; g.Fill(0.0); for (int i=0; i theta; theta.set_size(NPoints); vnl_vector phi; phi.set_size(NPoints); double C = sqrt(4*M_PI); phi(0) = 0.0; phi(NPoints-1) = 0.0; for(int i=0; i0 && i std::vector > QmitkFiberfoxView::MakeGradientList() { std::vector > retval; vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); // Add 0 vector for B0 int numB0 = ndirs/10; if (numB0==0) numB0=1; itk::Vector v; v.Fill(0.0); for (int i=0; i v; v[0] = U->get(0,i); v[1] = U->get(1,i); v[2] = U->get(2,i); retval.push_back(v); } return retval; } void QmitkFiberfoxView::OnAddBundle() { if (m_SelectedImage.IsNull()) return; mitk::DataStorage::SetOfObjects::ConstPointer children = GetDataStorage()->GetDerivations(m_SelectedImage); mitk::FiberBundleX::Pointer bundle = mitk::FiberBundleX::New(); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( bundle ); QString name = QString("Bundle_%1").arg(children->size()); node->SetName(name.toStdString()); m_SelectedBundles.push_back(node); UpdateGui(); GetDataStorage()->Add(node, m_SelectedImage); } void QmitkFiberfoxView::OnDrawROI() { if (m_SelectedBundles.empty()) OnAddBundle(); if (m_SelectedBundles.empty()) return; mitk::DataStorage::SetOfObjects::ConstPointer children = GetDataStorage()->GetDerivations(m_SelectedBundles.at(0)); mitk::PlanarEllipse::Pointer figure = mitk::PlanarEllipse::New(); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( figure ); node->SetBoolProperty("planarfigure.3drendering", true); QList nodes = this->GetDataManagerSelection(); for( int i=0; iSetSelected(false); m_SelectedFiducial = node; QString name = QString("Fiducial_%1").arg(children->size()); node->SetName(name.toStdString()); node->SetSelected(true); this->DisableCrosshairNavigation(); mitk::PlanarFigureInteractor::Pointer figureInteractor = dynamic_cast(node->GetDataInteractor().GetPointer()); if(figureInteractor.IsNull()) { figureInteractor = mitk::PlanarFigureInteractor::New(); us::Module* planarFigureModule = us::ModuleRegistry::GetModule( "MitkPlanarFigure" ); figureInteractor->LoadStateMachine("PlanarFigureInteraction.xml", planarFigureModule ); figureInteractor->SetEventConfig( "PlanarFigureConfig.xml", planarFigureModule ); figureInteractor->SetDataNode( node ); } UpdateGui(); GetDataStorage()->Add(node, m_SelectedBundles.at(0)); } bool CompareLayer(mitk::DataNode::Pointer i,mitk::DataNode::Pointer j) { int li = -1; i->GetPropertyValue("layer", li); int lj = -1; j->GetPropertyValue("layer", lj); return liGetSources(m_SelectedFiducial); for( mitk::DataStorage::SetOfObjects::const_iterator it = parents->begin(); it != parents->end(); ++it ) if(dynamic_cast((*it)->GetData())) m_SelectedBundles.push_back(*it); if (m_SelectedBundles.empty()) return; } FiberfoxParameters parameters = UpdateImageParameters(); for (unsigned int i=0; iGetDerivations(m_SelectedBundles.at(i)); std::vector< mitk::DataNode::Pointer > childVector; for( mitk::DataStorage::SetOfObjects::const_iterator it = children->begin(); it != children->end(); ++it ) childVector.push_back(*it); sort(childVector.begin(), childVector.end(), CompareLayer); vector< mitk::PlanarEllipse::Pointer > fib; vector< unsigned int > flip; float radius = 1; int count = 0; for( std::vector< mitk::DataNode::Pointer >::const_iterator it = childVector.begin(); it != childVector.end(); ++it ) { mitk::DataNode::Pointer node = *it; if ( node.IsNotNull() && dynamic_cast(node->GetData()) ) { mitk::PlanarEllipse* ellipse = dynamic_cast(node->GetData()); if (m_Controls->m_ConstantRadiusBox->isChecked()) { ellipse->SetTreatAsCircle(true); mitk::Point2D c = ellipse->GetControlPoint(0); mitk::Point2D p = ellipse->GetControlPoint(1); mitk::Vector2D v = p-c; if (count==0) { radius = v.GetVnlVector().magnitude(); ellipse->SetControlPoint(1, p); } else { v.Normalize(); v *= radius; ellipse->SetControlPoint(1, c+v); } } fib.push_back(ellipse); std::map::iterator it = m_DataNodeToPlanarFigureData.find(node.GetPointer()); if( it != m_DataNodeToPlanarFigureData.end() ) { QmitkPlanarFigureData& data = it->second; flip.push_back(data.m_Flipped); } else flip.push_back(0); } count++; } if (fib.size()>1) { parameters.m_FiberGen.m_Fiducials.push_back(fib); parameters.m_FiberGen.m_FlipList.push_back(flip); } else if (fib.size()>0) m_SelectedBundles.at(i)->SetData( mitk::FiberBundleX::New() ); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } itk::FibersFromPlanarFiguresFilter::Pointer filter = itk::FibersFromPlanarFiguresFilter::New(); filter->SetParameters(parameters.m_FiberGen); filter->Update(); vector< mitk::FiberBundleX::Pointer > fiberBundles = filter->GetFiberBundles(); for (unsigned int i=0; iSetData( fiberBundles.at(i) ); if (fiberBundles.at(i)->GetNumFibers()>50000) m_SelectedBundles.at(i)->SetVisibility(false); } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberfoxView::GenerateImage() { if (m_SelectedBundles.empty() && m_SelectedDWI.IsNull()) { mitk::Image::Pointer image = mitk::ImageGenerator::GenerateGradientImage( m_Controls->m_SizeX->value(), m_Controls->m_SizeY->value(), m_Controls->m_SizeZ->value(), m_Controls->m_SpacingX->value(), m_Controls->m_SpacingY->value(), m_Controls->m_SpacingZ->value()); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); node->SetName("Dummy"); unsigned int window = m_Controls->m_SizeX->value()*m_Controls->m_SizeY->value()*m_Controls->m_SizeZ->value(); unsigned int level = window/2; mitk::LevelWindow lw; lw.SetLevelWindow(level, window); node->SetProperty( "levelwindow", mitk::LevelWindowProperty::New( lw ) ); GetDataStorage()->Add(node); m_SelectedImage = node; mitk::BaseData::Pointer basedata = node->GetData(); if (basedata.IsNotNull()) { mitk::RenderingManager::GetInstance()->InitializeViews( basedata->GetTimeGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true ); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } UpdateGui(); } else if (!m_SelectedBundles.empty()) SimulateImageFromFibers(m_SelectedBundles.at(0)); else if (m_SelectedDWI.IsNotNull()) SimulateForExistingDwi(m_SelectedDWI); } void QmitkFiberfoxView::SimulateForExistingDwi(mitk::DataNode* imageNode) { if (!dynamic_cast*>(imageNode->GetData())) return; FiberfoxParameters parameters = UpdateImageParameters(); if (parameters.m_NoiseModel==NULL && parameters.m_SignalGen.m_Spikes==0 && parameters.m_SignalGen.m_FrequencyMap.IsNull() && parameters.m_SignalGen.m_KspaceLineOffset<=0.000001 && !parameters.m_SignalGen.m_DoAddGibbsRinging && !(parameters.m_SignalGen.m_EddyStrength>0) && parameters.m_SignalGen.m_CroppingFactor>0.999) { QMessageBox::information( NULL, "Simulation cancelled", "No valid artifact enabled! Motion artifacts and relaxation effects can NOT be added to an existing diffusion weighted image."); return; } mitk::DiffusionImage::Pointer diffImg = dynamic_cast*>(imageNode->GetData()); m_ArtifactsToDwiFilter = itk::AddArtifactsToDwiImageFilter< short >::New(); m_ArtifactsToDwiFilter->SetInput(diffImg->GetVectorImage()); parameters.m_Misc.m_ParentNode = imageNode; m_ArtifactsToDwiFilter->SetParameters(parameters); m_Worker.m_FilterType = 1; m_Thread.start(QThread::LowestPriority); } void QmitkFiberfoxView::SimulateImageFromFibers(mitk::DataNode* fiberNode) { mitk::FiberBundleX::Pointer fiberBundle = dynamic_cast(fiberNode->GetData()); if (fiberBundle->GetNumFibers()<=0) return; FiberfoxParameters parameters = UpdateImageParameters(); m_TractsToDwiFilter = itk::TractsToDWIImageFilter< short >::New(); parameters.m_Misc.m_ParentNode = fiberNode; if (m_SelectedDWI.IsNotNull()) { mitk::DiffusionImage::Pointer diffImg = dynamic_cast*>(m_SelectedDWI->GetData()); bool doSampling = false; + bool ok = true; for (unsigned int i=0; i* >(parameters.m_FiberModelList[i]) ) + { doSampling = true; + RawShModel* model = dynamic_cast< RawShModel* >(parameters.m_FiberModelList[i]); + if (ok==true) + ok = model->SampleKernels(diffImg, parameters.m_SignalGen.m_MaskImage); + } for (unsigned int i=0; i* >(parameters.m_NonFiberModelList[i]) ) - doSampling = true; - - // sample prototype signals - if ( doSampling ) - { - const int shOrder = 2; - - typedef itk::DiffusionTensor3DReconstructionImageFilter< short, short, double > TensorReconstructionImageFilterType; - TensorReconstructionImageFilterType::Pointer filter = TensorReconstructionImageFilterType::New(); - filter->SetGradientImage( diffImg->GetDirections(), diffImg->GetVectorImage() ); - filter->SetBValue(diffImg->GetReferenceBValue()); - filter->Update(); - itk::Image< itk::DiffusionTensor3D< double >, 3 >::Pointer tensorImage = filter->GetOutput(); - - const int NumCoeffs = (shOrder*shOrder + shOrder + 2)/2 + shOrder; - typedef itk::AnalyticalDiffusionQballReconstructionImageFilter QballFilterType; - QballFilterType::Pointer qballfilter = QballFilterType::New(); - qballfilter->SetGradientImage( diffImg->GetDirections(), diffImg->GetVectorImage() ); - qballfilter->SetBValue(diffImg->GetReferenceBValue()); - qballfilter->SetLambda(0.006); - qballfilter->SetNormalizationMethod(QballFilterType::QBAR_RAW_SIGNAL); - qballfilter->Update(); - QballFilterType::CoefficientImageType::Pointer itkFeatureImage = qballfilter->GetCoefficientImage(); - - itk::AdcImageFilter< short, double >::Pointer adcFilter = itk::AdcImageFilter< short, double >::New(); - adcFilter->SetInput(diffImg->GetVectorImage()); - adcFilter->SetGradientDirections(diffImg->GetDirections()); - adcFilter->SetB_value(diffImg->GetReferenceBValue()); - adcFilter->Update(); - ItkDoubleImgType::Pointer adcImage = adcFilter->GetOutput(); - - int b0Index; - for (unsigned int i=0; iGetDirectionsWithoutMeasurementFrame()->Size(); i++) - if ( diffImg->GetDirectionsWithoutMeasurementFrame()->GetElement(i).magnitude()<0.001 ) - { - b0Index = i; - break; - } - - double max = 0; - { - itk::ImageRegionIterator< itk::VectorImage< short, 3 > > it(diffImg->GetVectorImage(), diffImg->GetVectorImage()->GetLargestPossibleRegion()); - while(!it.IsAtEnd()) - { - if (parameters.m_SignalGen.m_MaskImage.IsNotNull() && parameters.m_SignalGen.m_MaskImage->GetPixel(it.GetIndex())<=0) - { - ++it; - continue; - } - if (it.Get()[b0Index]>max) - max = it.Get()[b0Index]; - ++it; - } - } - - MITK_INFO << "Sampling signal kernels."; - - itk::ImageRegionIterator< itk::Image< itk::DiffusionTensor3D< double >, 3 > > it(tensorImage, tensorImage->GetLargestPossibleRegion()); - while(!it.IsAtEnd()) - { - bool skipPixel = false; - if (parameters.m_SignalGen.m_MaskImage.IsNotNull() && parameters.m_SignalGen.m_MaskImage->GetPixel(it.GetIndex())<=0) - { - ++it; - continue; - } - for (unsigned int i=0; iGetDirections()->Size(); i++) - { - if (diffImg->GetVectorImage()->GetPixel(it.GetIndex())[i]!=diffImg->GetVectorImage()->GetPixel(it.GetIndex())[i] || diffImg->GetVectorImage()->GetPixel(it.GetIndex())[i]<=0 || diffImg->GetVectorImage()->GetPixel(it.GetIndex())[i]>diffImg->GetVectorImage()->GetPixel(it.GetIndex())[b0Index]) - { - skipPixel = true; - break; - } - } - if (skipPixel) - { - ++it; - continue; - } - - typedef itk::DiffusionTensor3D TensorType; - TensorType::EigenValuesArrayType eigenvalues; - TensorType::EigenVectorsMatrixType eigenvectors; - TensorType tensor = it.Get(); - double FA = tensor.GetFractionalAnisotropy(); - double ADC = adcImage->GetPixel(it.GetIndex()); - QballFilterType::CoefficientImageType::PixelType itkv = itkFeatureImage->GetPixel(it.GetIndex()); - vnl_vector_fixed< double, NumCoeffs > coeffs; - for (unsigned int c=0; c* model = dynamic_cast< RawShModel* >(parameters.m_FiberModelList[i]); - if ( model && model->GetMaxNumKernels()>model->GetNumberOfKernels() && FA>model->GetFaRange().first && FAGetFaRange().second && ADC>model->GetAdcRange().first && ADCGetAdcRange().second) - { - if (model->SetShCoefficients( coeffs, (double)diffImg->GetVectorImage()->GetPixel(it.GetIndex())[b0Index]/max )) - { - tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); - itk::Vector dir; - dir[0] = eigenvectors(2, 0); - dir[1] = eigenvectors(2, 1); - dir[2] = eigenvectors(2, 2); - model->m_PrototypeMaxDirection.push_back(dir); - MITK_INFO << "WM KERNEL: " << it.GetIndex() << " (" << model->GetNumberOfKernels() << ")"; - } - } - } - for (unsigned int i=0; i* model = dynamic_cast< RawShModel* >(parameters.m_NonFiberModelList[i]); - if ( model && model->GetMaxNumKernels()>model->GetNumberOfKernels() && FA>model->GetFaRange().first && FAGetFaRange().second && ADC>model->GetAdcRange().first && ADCGetAdcRange().second) - { - if (model->SetShCoefficients( coeffs, (double)diffImg->GetVectorImage()->GetPixel(it.GetIndex())[b0Index]/max )) - MITK_INFO << "CSF/GM KERNEL: " << it.GetIndex() << " (" << model->GetNumberOfKernels() << ")"; - } - } - ++it; - } - - for (unsigned int i=0; i* model = dynamic_cast< RawShModel* >(parameters.m_FiberModelList[i]); - if ( model && model->GetNumberOfKernels()<=0 ) - { - QMessageBox::information( NULL, "Simulation cancelled", "No suitable voxels found for fiber compartment "+QString::number(i)); - return; - } - } - for (unsigned int i=0; i* model = dynamic_cast< RawShModel* >(parameters.m_NonFiberModelList[i]); - if ( model && model->GetNumberOfKernels()<=0 ) - { - QMessageBox::information( NULL, "Simulation cancelled", "No suitable voxels found for non-fiber compartment "+QString::number(i)); - return; - } + if (ok==true) + ok = model->SampleKernels(diffImg, parameters.m_SignalGen.m_MaskImage); } + +// // sample prototype signals +// if ( doSampling ) +// { +// const int shOrder = 2; + +// typedef itk::DiffusionTensor3DReconstructionImageFilter< short, short, double > TensorReconstructionImageFilterType; +// TensorReconstructionImageFilterType::Pointer filter = TensorReconstructionImageFilterType::New(); +// filter->SetGradientImage( diffImg->GetDirections(), diffImg->GetVectorImage() ); +// filter->SetBValue(diffImg->GetReferenceBValue()); +// filter->Update(); +// itk::Image< itk::DiffusionTensor3D< double >, 3 >::Pointer tensorImage = filter->GetOutput(); + +// const int NumCoeffs = (shOrder*shOrder + shOrder + 2)/2 + shOrder; +// typedef itk::AnalyticalDiffusionQballReconstructionImageFilter QballFilterType; +// QballFilterType::Pointer qballfilter = QballFilterType::New(); +// qballfilter->SetGradientImage( diffImg->GetDirections(), diffImg->GetVectorImage() ); +// qballfilter->SetBValue(diffImg->GetReferenceBValue()); +// qballfilter->SetLambda(0.006); +// qballfilter->SetNormalizationMethod(QballFilterType::QBAR_RAW_SIGNAL); +// qballfilter->Update(); +// QballFilterType::CoefficientImageType::Pointer itkFeatureImage = qballfilter->GetCoefficientImage(); + +// itk::AdcImageFilter< short, double >::Pointer adcFilter = itk::AdcImageFilter< short, double >::New(); +// adcFilter->SetInput(diffImg->GetVectorImage()); +// adcFilter->SetGradientDirections(diffImg->GetDirections()); +// adcFilter->SetB_value(diffImg->GetReferenceBValue()); +// adcFilter->Update(); +// ItkDoubleImgType::Pointer adcImage = adcFilter->GetOutput(); +// } + if (!ok) + { + QMessageBox::information( NULL, "Simulation cancelled", "No valid prototype signals could be sampled."); + return; } } else if ( m_Controls->m_Compartment1Box->currentIndex()==3 || m_Controls->m_Compartment3Box->currentIndex()==3 || m_Controls->m_Compartment4Box->currentIndex()==4 ) { QMessageBox::information( NULL, "Simulation cancelled", "Prototype signal but no diffusion-weighted image selected to sample signal from."); return; } m_TractsToDwiFilter->SetParameters(parameters); m_TractsToDwiFilter->SetFiberBundle(fiberBundle); m_Worker.m_FilterType = 0; m_Thread.start(QThread::LowestPriority); } void QmitkFiberfoxView::ApplyTransform() { vector< mitk::DataNode::Pointer > selectedBundles; for(unsigned int i=0; iGetDerivations(m_SelectedImages.at(i)); for( mitk::DataStorage::SetOfObjects::const_iterator it = derivations->begin(); it != derivations->end(); ++it ) { mitk::DataNode::Pointer fibNode = *it; if ( fibNode.IsNotNull() && dynamic_cast(fibNode->GetData()) ) selectedBundles.push_back(fibNode); } } if (selectedBundles.empty()) selectedBundles = m_SelectedBundles2; if (!selectedBundles.empty()) { for (std::vector::const_iterator it = selectedBundles.begin(); it!=selectedBundles.end(); ++it) { mitk::FiberBundleX::Pointer fib = dynamic_cast((*it)->GetData()); fib->RotateAroundAxis(m_Controls->m_XrotBox->value(), m_Controls->m_YrotBox->value(), m_Controls->m_ZrotBox->value()); fib->TranslateFibers(m_Controls->m_XtransBox->value(), m_Controls->m_YtransBox->value(), m_Controls->m_ZtransBox->value()); fib->ScaleFibers(m_Controls->m_XscaleBox->value(), m_Controls->m_YscaleBox->value(), m_Controls->m_ZscaleBox->value()); // handle child fiducials if (m_Controls->m_IncludeFiducials->isChecked()) { mitk::DataStorage::SetOfObjects::ConstPointer derivations = GetDataStorage()->GetDerivations(*it); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = derivations->begin(); it2 != derivations->end(); ++it2 ) { mitk::DataNode::Pointer fiducialNode = *it2; if ( fiducialNode.IsNotNull() && dynamic_cast(fiducialNode->GetData()) ) { mitk::PlanarEllipse* pe = dynamic_cast(fiducialNode->GetData()); mitk::BaseGeometry* geom = pe->GetGeometry(); // translate mitk::Vector3D world; world[0] = m_Controls->m_XtransBox->value(); world[1] = m_Controls->m_YtransBox->value(); world[2] = m_Controls->m_ZtransBox->value(); geom->Translate(world); // calculate rotation matrix double x = m_Controls->m_XrotBox->value()*M_PI/180; double y = m_Controls->m_YrotBox->value()*M_PI/180; double z = m_Controls->m_ZrotBox->value()*M_PI/180; itk::Matrix< double, 3, 3 > rotX; rotX.SetIdentity(); rotX[1][1] = cos(x); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(x); rotX[2][1] = -rotX[1][2]; itk::Matrix< double, 3, 3 > rotY; rotY.SetIdentity(); rotY[0][0] = cos(y); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(y); rotY[2][0] = -rotY[0][2]; itk::Matrix< double, 3, 3 > rotZ; rotZ.SetIdentity(); rotZ[0][0] = cos(z); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(z); rotZ[1][0] = -rotZ[0][1]; itk::Matrix< double, 3, 3 > rot = rotZ*rotY*rotX; // transform control point coordinate into geometry translation geom->SetOrigin(pe->GetWorldControlPoint(0)); mitk::Point2D cp; cp.Fill(0.0); pe->SetControlPoint(0, cp); // rotate fiducial geom->GetIndexToWorldTransform()->SetMatrix(rot*geom->GetIndexToWorldTransform()->GetMatrix()); // implicit translation mitk::Vector3D trans; trans[0] = geom->GetOrigin()[0]-fib->GetGeometry()->GetCenter()[0]; trans[1] = geom->GetOrigin()[1]-fib->GetGeometry()->GetCenter()[1]; trans[2] = geom->GetOrigin()[2]-fib->GetGeometry()->GetCenter()[2]; mitk::Vector3D newWc = rot*trans; newWc = newWc-trans; geom->Translate(newWc); pe->Modified(); } } } } } else { for (unsigned int i=0; i(m_SelectedFiducials.at(i)->GetData()); mitk::BaseGeometry* geom = pe->GetGeometry(); // translate mitk::Vector3D world; world[0] = m_Controls->m_XtransBox->value(); world[1] = m_Controls->m_YtransBox->value(); world[2] = m_Controls->m_ZtransBox->value(); geom->Translate(world); // calculate rotation matrix double x = m_Controls->m_XrotBox->value()*M_PI/180; double y = m_Controls->m_YrotBox->value()*M_PI/180; double z = m_Controls->m_ZrotBox->value()*M_PI/180; itk::Matrix< double, 3, 3 > rotX; rotX.SetIdentity(); rotX[1][1] = cos(x); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(x); rotX[2][1] = -rotX[1][2]; itk::Matrix< double, 3, 3 > rotY; rotY.SetIdentity(); rotY[0][0] = cos(y); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(y); rotY[2][0] = -rotY[0][2]; itk::Matrix< double, 3, 3 > rotZ; rotZ.SetIdentity(); rotZ[0][0] = cos(z); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(z); rotZ[1][0] = -rotZ[0][1]; itk::Matrix< double, 3, 3 > rot = rotZ*rotY*rotX; // transform control point coordinate into geometry translation geom->SetOrigin(pe->GetWorldControlPoint(0)); mitk::Point2D cp; cp.Fill(0.0); pe->SetControlPoint(0, cp); // rotate fiducial geom->GetIndexToWorldTransform()->SetMatrix(rot*geom->GetIndexToWorldTransform()->GetMatrix()); pe->Modified(); } if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberfoxView::CopyBundles() { if ( m_SelectedBundles.size()<1 ){ QMessageBox::information( NULL, "Warning", "Select at least one fiber bundle!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least one fiber bundle!"; return; } for (std::vector::const_iterator it = m_SelectedBundles.begin(); it!=m_SelectedBundles.end(); ++it) { // find parent image mitk::DataNode::Pointer parentNode; mitk::DataStorage::SetOfObjects::ConstPointer parentImgs = GetDataStorage()->GetSources(*it); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = parentImgs->begin(); it2 != parentImgs->end(); ++it2 ) { mitk::DataNode::Pointer pImgNode = *it2; if ( pImgNode.IsNotNull() && dynamic_cast(pImgNode->GetData()) ) { parentNode = pImgNode; break; } } mitk::FiberBundleX::Pointer fib = dynamic_cast((*it)->GetData()); mitk::FiberBundleX::Pointer newBundle = fib->GetDeepCopy(); QString name((*it)->GetName().c_str()); name += "_copy"; mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); if (parentNode.IsNotNull()) GetDataStorage()->Add(fbNode, parentNode); else GetDataStorage()->Add(fbNode); // copy child fiducials if (m_Controls->m_IncludeFiducials->isChecked()) { mitk::DataStorage::SetOfObjects::ConstPointer derivations = GetDataStorage()->GetDerivations(*it); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = derivations->begin(); it2 != derivations->end(); ++it2 ) { mitk::DataNode::Pointer fiducialNode = *it2; if ( fiducialNode.IsNotNull() && dynamic_cast(fiducialNode->GetData()) ) { mitk::PlanarEllipse::Pointer pe = mitk::PlanarEllipse::New(); pe->DeepCopy(dynamic_cast(fiducialNode->GetData())); mitk::DataNode::Pointer newNode = mitk::DataNode::New(); newNode->SetData(pe); newNode->SetName(fiducialNode->GetName()); newNode->SetBoolProperty("planarfigure.3drendering", true); GetDataStorage()->Add(newNode, fbNode); } } } } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberfoxView::JoinBundles() { if ( m_SelectedBundles.size()<2 ){ QMessageBox::information( NULL, "Warning", "Select at least two fiber bundles!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!"; return; } std::vector::const_iterator it = m_SelectedBundles.begin(); mitk::FiberBundleX::Pointer newBundle = dynamic_cast((*it)->GetData()); QString name(""); name += QString((*it)->GetName().c_str()); ++it; for (; it!=m_SelectedBundles.end(); ++it) { newBundle = newBundle->AddBundle(dynamic_cast((*it)->GetData())); name += "+"+QString((*it)->GetName().c_str()); } mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberfoxView::UpdateGui() { m_Controls->m_FiberBundleLabel->setText("mandatory"); m_Controls->m_GeometryFrame->setEnabled(true); m_Controls->m_GeometryMessage->setVisible(false); m_Controls->m_DiffusionPropsMessage->setVisible(false); m_Controls->m_FiberGenMessage->setVisible(true); m_Controls->m_TransformBundlesButton->setEnabled(false); m_Controls->m_CopyBundlesButton->setEnabled(false); m_Controls->m_GenerateFibersButton->setEnabled(false); m_Controls->m_FlipButton->setEnabled(false); m_Controls->m_CircleButton->setEnabled(false); m_Controls->m_BvalueBox->setEnabled(true); m_Controls->m_NumGradientsBox->setEnabled(true); m_Controls->m_JoinBundlesButton->setEnabled(false); m_Controls->m_AlignOnGrid->setEnabled(false); if (m_SelectedFiducial.IsNotNull()) { m_Controls->m_TransformBundlesButton->setEnabled(true); m_Controls->m_FlipButton->setEnabled(true); m_Controls->m_AlignOnGrid->setEnabled(true); } if (m_SelectedImage.IsNotNull() || !m_SelectedBundles.empty()) { m_Controls->m_TransformBundlesButton->setEnabled(true); m_Controls->m_CircleButton->setEnabled(true); m_Controls->m_FiberGenMessage->setVisible(false); m_Controls->m_AlignOnGrid->setEnabled(true); } if (m_MaskImageNode.IsNotNull() || m_SelectedImage.IsNotNull()) { m_Controls->m_GeometryMessage->setVisible(true); m_Controls->m_GeometryFrame->setEnabled(false); } if (m_SelectedDWI.IsNotNull()) { m_Controls->m_DiffusionPropsMessage->setVisible(true); m_Controls->m_BvalueBox->setEnabled(false); m_Controls->m_NumGradientsBox->setEnabled(false); m_Controls->m_GeometryMessage->setVisible(true); m_Controls->m_GeometryFrame->setEnabled(false); } if (!m_SelectedBundles.empty()) { m_Controls->m_CopyBundlesButton->setEnabled(true); m_Controls->m_GenerateFibersButton->setEnabled(true); m_Controls->m_FiberBundleLabel->setText(m_SelectedBundles.at(0)->GetName().c_str()); if (m_SelectedBundles.size()>1) m_Controls->m_JoinBundlesButton->setEnabled(true); } } void QmitkFiberfoxView::OnSelectionChanged( berry::IWorkbenchPart::Pointer, const QList& nodes ) { m_SelectedBundles2.clear(); m_SelectedImages.clear(); m_SelectedFiducials.clear(); m_SelectedFiducial = NULL; m_SelectedBundles.clear(); m_SelectedImage = NULL; m_SelectedDWI = NULL; m_MaskImageNode = NULL; m_Controls->m_TissueMaskLabel->setText("optional"); // iterate all selected objects, adjust warning visibility for( int i=0; i*>(node->GetData()) ) { m_SelectedDWI = node; m_SelectedImage = node; m_SelectedImages.push_back(node); } else if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { m_SelectedImages.push_back(node); m_SelectedImage = node; bool isbinary = false; node->GetPropertyValue("binary", isbinary); if (isbinary) { m_MaskImageNode = node; m_Controls->m_TissueMaskLabel->setText(m_MaskImageNode->GetName().c_str()); } } else if ( node.IsNotNull() && dynamic_cast(node->GetData()) ) { m_SelectedBundles2.push_back(node); if (m_Controls->m_RealTimeFibers->isChecked()) { m_SelectedBundles.push_back(node); mitk::FiberBundleX::Pointer newFib = dynamic_cast(node->GetData()); if (newFib->GetNumFibers()!=m_Controls->m_FiberDensityBox->value()) GenerateFibers(); } else m_SelectedBundles.push_back(node); } else if ( node.IsNotNull() && dynamic_cast(node->GetData()) ) { m_SelectedFiducials.push_back(node); m_SelectedFiducial = node; m_SelectedBundles.clear(); mitk::DataStorage::SetOfObjects::ConstPointer parents = GetDataStorage()->GetSources(node); for( mitk::DataStorage::SetOfObjects::const_iterator it = parents->begin(); it != parents->end(); ++it ) { mitk::DataNode::Pointer pNode = *it; if ( pNode.IsNotNull() && dynamic_cast(pNode->GetData()) ) m_SelectedBundles.push_back(pNode); } } } UpdateGui(); } void QmitkFiberfoxView::EnableCrosshairNavigation() { MITK_DEBUG << "EnableCrosshairNavigation"; // enable the crosshair navigation if (mitk::ILinkedRenderWindowPart* linkedRenderWindow = dynamic_cast(this->GetRenderWindowPart())) { MITK_DEBUG << "enabling linked navigation"; linkedRenderWindow->EnableLinkedNavigation(true); // linkedRenderWindow->EnableSlicingPlanes(true); } if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::DisableCrosshairNavigation() { MITK_DEBUG << "DisableCrosshairNavigation"; // disable the crosshair navigation during the drawing if (mitk::ILinkedRenderWindowPart* linkedRenderWindow = dynamic_cast(this->GetRenderWindowPart())) { MITK_DEBUG << "disabling linked navigation"; linkedRenderWindow->EnableLinkedNavigation(false); // linkedRenderWindow->EnableSlicingPlanes(false); } } void QmitkFiberfoxView::NodeRemoved(const mitk::DataNode* node) { mitk::DataNode* nonConstNode = const_cast(node); std::map::iterator it = m_DataNodeToPlanarFigureData.find(nonConstNode); if (dynamic_cast(node->GetData())) { m_SelectedBundles.clear(); m_SelectedBundles2.clear(); } else if (dynamic_cast(node->GetData())) m_SelectedImages.clear(); if( it != m_DataNodeToPlanarFigureData.end() ) { QmitkPlanarFigureData& data = it->second; // remove observers data.m_Figure->RemoveObserver( data.m_EndPlacementObserverTag ); data.m_Figure->RemoveObserver( data.m_SelectObserverTag ); data.m_Figure->RemoveObserver( data.m_StartInteractionObserverTag ); data.m_Figure->RemoveObserver( data.m_EndInteractionObserverTag ); m_DataNodeToPlanarFigureData.erase( it ); } } void QmitkFiberfoxView::NodeAdded( const mitk::DataNode* node ) { // add observer for selection in renderwindow mitk::PlanarFigure* figure = dynamic_cast(node->GetData()); bool isPositionMarker (false); node->GetBoolProperty("isContourMarker", isPositionMarker); if( figure && !isPositionMarker ) { MITK_DEBUG << "figure added. will add interactor if needed."; mitk::PlanarFigureInteractor::Pointer figureInteractor = dynamic_cast(node->GetDataInteractor().GetPointer()); mitk::DataNode* nonConstNode = const_cast( node ); if(figureInteractor.IsNull()) { figureInteractor = mitk::PlanarFigureInteractor::New(); us::Module* planarFigureModule = us::ModuleRegistry::GetModule( "MitkPlanarFigure" ); figureInteractor->LoadStateMachine("PlanarFigureInteraction.xml", planarFigureModule ); figureInteractor->SetEventConfig( "PlanarFigureConfig.xml", planarFigureModule ); figureInteractor->SetDataNode( nonConstNode ); } MITK_DEBUG << "will now add observers for planarfigure"; QmitkPlanarFigureData data; data.m_Figure = figure; // // add observer for event when figure has been placed typedef itk::SimpleMemberCommand< QmitkFiberfoxView > SimpleCommandType; // SimpleCommandType::Pointer initializationCommand = SimpleCommandType::New(); // initializationCommand->SetCallbackFunction( this, &QmitkFiberfoxView::PlanarFigureInitialized ); // data.m_EndPlacementObserverTag = figure->AddObserver( mitk::EndPlacementPlanarFigureEvent(), initializationCommand ); // add observer for event when figure is picked (selected) typedef itk::MemberCommand< QmitkFiberfoxView > MemberCommandType; MemberCommandType::Pointer selectCommand = MemberCommandType::New(); selectCommand->SetCallbackFunction( this, &QmitkFiberfoxView::PlanarFigureSelected ); data.m_SelectObserverTag = figure->AddObserver( mitk::SelectPlanarFigureEvent(), selectCommand ); // add observer for event when interaction with figure starts SimpleCommandType::Pointer startInteractionCommand = SimpleCommandType::New(); startInteractionCommand->SetCallbackFunction( this, &QmitkFiberfoxView::DisableCrosshairNavigation); data.m_StartInteractionObserverTag = figure->AddObserver( mitk::StartInteractionPlanarFigureEvent(), startInteractionCommand ); // add observer for event when interaction with figure starts SimpleCommandType::Pointer endInteractionCommand = SimpleCommandType::New(); endInteractionCommand->SetCallbackFunction( this, &QmitkFiberfoxView::EnableCrosshairNavigation); data.m_EndInteractionObserverTag = figure->AddObserver( mitk::EndInteractionPlanarFigureEvent(), endInteractionCommand ); m_DataNodeToPlanarFigureData[nonConstNode] = data; } } void QmitkFiberfoxView::PlanarFigureSelected( itk::Object* object, const itk::EventObject& ) { mitk::TNodePredicateDataType::Pointer isPf = mitk::TNodePredicateDataType::New(); mitk::DataStorage::SetOfObjects::ConstPointer allPfs = this->GetDataStorage()->GetSubset( isPf ); for ( mitk::DataStorage::SetOfObjects::const_iterator it = allPfs->begin(); it!=allPfs->end(); ++it) { mitk::DataNode* node = *it; if( node->GetData() == object ) { node->SetSelected(true); m_SelectedFiducial = node; } else node->SetSelected(false); } UpdateGui(); this->RequestRenderWindowUpdate(); } void QmitkFiberfoxView::SetFocus() { m_Controls->m_CircleButton->setFocus(); } void QmitkFiberfoxView::SetOutputPath() { // SELECT FOLDER DIALOG string outputPath = QFileDialog::getExistingDirectory(NULL, "Save images to...", QString(outputPath.c_str())).toStdString(); if (outputPath.empty()) m_Controls->m_SavePathEdit->setText("-"); else { outputPath += "/"; m_Controls->m_SavePathEdit->setText(QString(outputPath.c_str())); } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.h b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.h index dbe58f859b..be24a54ab1 100755 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.h @@ -1,209 +1,210 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include "ui_QmitkFiberfoxViewControls.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include +#include #include /*! \brief View for fiber based diffusion software phantoms (Fiberfox). See "Fiberfox: Facilitating the creation of realistic white matter software phantoms" (DOI: 10.1002/mrm.25045) for details. \sa QmitkFunctionality \ingroup Functionalities */ // Forward Qt class declarations using namespace std; class QmitkFiberfoxView; class QmitkFiberfoxWorker : public QObject { Q_OBJECT public: QmitkFiberfoxWorker(QmitkFiberfoxView* view); int m_FilterType; public slots: void run(); private: QmitkFiberfoxView* m_View; }; class QmitkFiberfoxView : public QmitkAbstractView { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: static const string VIEW_ID; QmitkFiberfoxView(); virtual ~QmitkFiberfoxView(); virtual void CreateQtPartControl(QWidget *parent); void SetFocus(); typedef itk::Image ItkDoubleImgType; typedef itk::Image ItkUcharImgType; typedef itk::Vector GradientType; typedef vector GradientListType; template vector > MakeGradientList(); protected slots: void SetOutputPath(); ///< path where image is automatically saved to after the simulation is finished void LoadParameters(); ///< load fiberfox parameters void SaveParameters(); ///< save fiberfox parameters void BeforeThread(); void AfterThread(); void KillThread(); ///< abort simulation void UpdateSimulationStatus(); ///< print simulation progress and satus messages void OnDrawROI(); ///< adds new ROI, handles interactors etc. void OnAddBundle(); ///< adds new fiber bundle to datastorage void OnFlipButton(); ///< negate one coordinate of the fiber waypoints in the selcted planar figure. needed in case of unresolvable twists void GenerateFibers(); ///< generate fibers from the selected ROIs void GenerateImage(); ///< start image simulation void JoinBundles(); ///< merges selcted fiber bundles into one void CopyBundles(); ///< add copy of the selected bundle to the datamanager void ApplyTransform(); ///< rotate and shift selected bundles void AlignOnGrid(); ///< shift selected fiducials to nearest voxel center void Comp1ModelFrameVisibility(int index); ///< only show parameters of selected signal model for compartment 1 void Comp2ModelFrameVisibility(int index); ///< only show parameters of selected signal model for compartment 2 void Comp3ModelFrameVisibility(int index); ///< only show parameters of selected signal model for compartment 3 void Comp4ModelFrameVisibility(int index); ///< only show parameters of selected signal model for compartment 4 void ShowAdvancedOptions(int state); /** update fibers if any parameter changes */ void OnFiberDensityChanged(int value); void OnFiberSamplingChanged(double value); void OnTensionChanged(double value); void OnContinuityChanged(double value); void OnBiasChanged(double value); void OnVarianceChanged(double value); void OnDistributionChanged(int value); void OnConstantRadius(int value); /** update GUI elements */ void OnAddNoise(int value); void OnAddGhosts(int value); void OnAddDistortions(int value); void OnAddEddy(int value); void OnAddSpikes(int value); void OnAddAliasing(int value); void OnAddMotion(int value); protected: /// \brief called by QmitkFunctionality when DataManager's selection has changed virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer, const QList&); GradientListType GenerateHalfShell(int NPoints); ///< generate vectors distributed over the halfsphere Ui::QmitkFiberfoxViewControls* m_Controls; void SimulateForExistingDwi(mitk::DataNode* imageNode); ///< add artifacts to existing diffusion weighted image void SimulateImageFromFibers(mitk::DataNode* fiberNode); ///< simulate new diffusion weighted image template< class ScalarType > FiberfoxParameters< ScalarType > UpdateImageParameters(); ///< update fiberfox paramater object (template parameter defines noise model type) void UpdateGui(); ///< enable/disbale buttons etc. according to current datamanager selection void PlanarFigureSelected( itk::Object* object, const itk::EventObject& ); void EnableCrosshairNavigation(); ///< enable crosshair navigation if planar figure interaction ends void DisableCrosshairNavigation(); ///< disable crosshair navigation if planar figure interaction starts void NodeAdded( const mitk::DataNode* node ); ///< add observers void NodeRemoved(const mitk::DataNode* node); ///< remove observers /** structure to keep track of planar figures and observers */ struct QmitkPlanarFigureData { QmitkPlanarFigureData() : m_Figure(0) , m_EndPlacementObserverTag(0) , m_SelectObserverTag(0) , m_StartInteractionObserverTag(0) , m_EndInteractionObserverTag(0) , m_Flipped(0) { } mitk::PlanarFigure* m_Figure; unsigned int m_EndPlacementObserverTag; unsigned int m_SelectObserverTag; unsigned int m_StartInteractionObserverTag; unsigned int m_EndInteractionObserverTag; unsigned int m_Flipped; }; std::map m_DataNodeToPlanarFigureData; ///< map each planar figure uniquely to a QmitkPlanarFigureData mitk::DataNode::Pointer m_SelectedFiducial; ///< selected planar ellipse mitk::DataNode::Pointer m_SelectedImage; mitk::DataNode::Pointer m_SelectedDWI; vector< mitk::DataNode::Pointer > m_SelectedBundles; vector< mitk::DataNode::Pointer > m_SelectedBundles2; vector< mitk::DataNode::Pointer > m_SelectedFiducials; vector< mitk::DataNode::Pointer > m_SelectedImages; mitk::DataNode::Pointer m_MaskImageNode; QString m_ParameterFile; ///< parameter file name // GUI thread QmitkFiberfoxWorker m_Worker; ///< runs filter QThread m_Thread; ///< worker thread bool m_ThreadIsRunning; QTimer* m_SimulationTimer; QTime m_SimulationTime; QString m_SimulationStatusText; /** Image filters that do all the simulations. */ itk::TractsToDWIImageFilter< short >::Pointer m_TractsToDwiFilter; itk::AddArtifactsToDwiImageFilter< short >::Pointer m_ArtifactsToDwiFilter; friend class QmitkFiberfoxWorker; };