diff --git a/Modules/Core/src/DataManagement/mitkDataStorage.cpp b/Modules/Core/src/DataManagement/mitkDataStorage.cpp index b3f33ce251..0b18cd2c6b 100644 --- a/Modules/Core/src/DataManagement/mitkDataStorage.cpp +++ b/Modules/Core/src/DataManagement/mitkDataStorage.cpp @@ -1,506 +1,518 @@ /*=================================================================== 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 "mitkDataStorage.h" #include "itkCommand.h" #include "itkMutexLockHolder.h" #include "mitkDataNode.h" #include "mitkGroupTagProperty.h" #include "mitkImage.h" #include "mitkNodePredicateBase.h" #include "mitkNodePredicateProperty.h" #include "mitkProperties.h" mitk::DataStorage::DataStorage() : itk::Object(), m_BlockNodeModifiedEvents(false) { } mitk::DataStorage::~DataStorage() { ///// we can not call GetAll() in destructor, because it is implemented in a subclass // SetOfObjects::ConstPointer all = this->GetAll(); // for (SetOfObjects::ConstIterator it = all->Begin(); it != all->End(); ++it) // this->RemoveListeners(it->Value()); // m_NodeModifiedObserverTags.clear(); // m_NodeDeleteObserverTags.clear(); } void mitk::DataStorage::Add(mitk::DataNode *node, mitk::DataNode *parent) { mitk::DataStorage::SetOfObjects::Pointer parents = mitk::DataStorage::SetOfObjects::New(); if (parent != NULL) //< Return empty set if parent is null parents->InsertElement(0, parent); this->Add(node, parents); } void mitk::DataStorage::Remove(const mitk::DataStorage::SetOfObjects *nodes) { if (nodes == NULL) return; for (mitk::DataStorage::SetOfObjects::ConstIterator it = nodes->Begin(); it != nodes->End(); it++) this->Remove(it.Value()); } mitk::DataStorage::SetOfObjects::ConstPointer mitk::DataStorage::GetSubset(const NodePredicateBase *condition) const { mitk::DataStorage::SetOfObjects::ConstPointer result = this->FilterSetOfObjects(this->GetAll(), condition); return result; } mitk::DataNode *mitk::DataStorage::GetNamedNode(const char *name) const { if (name == NULL) return NULL; mitk::StringProperty::Pointer s(mitk::StringProperty::New(name)); mitk::NodePredicateProperty::Pointer p = mitk::NodePredicateProperty::New("name", s); mitk::DataStorage::SetOfObjects::ConstPointer rs = this->GetSubset(p); if (rs->Size() >= 1) return rs->GetElement(0); else return NULL; } mitk::DataNode *mitk::DataStorage::GetNode(const NodePredicateBase *condition) const { if (condition == NULL) return NULL; mitk::DataStorage::SetOfObjects::ConstPointer rs = this->GetSubset(condition); if (rs->Size() >= 1) return rs->GetElement(0); else return NULL; } mitk::DataNode *mitk::DataStorage::GetNamedDerivedNode(const char *name, const mitk::DataNode *sourceNode, bool onlyDirectDerivations) const { if (name == NULL) return NULL; mitk::StringProperty::Pointer s(mitk::StringProperty::New(name)); mitk::NodePredicateProperty::Pointer p = mitk::NodePredicateProperty::New("name", s); mitk::DataStorage::SetOfObjects::ConstPointer rs = this->GetDerivations(sourceNode, p, onlyDirectDerivations); if (rs->Size() >= 1) return rs->GetElement(0); else return NULL; } void mitk::DataStorage::PrintSelf(std::ostream &os, itk::Indent indent) const { // Superclass::PrintSelf(os, indent); mitk::DataStorage::SetOfObjects::ConstPointer all = this->GetAll(); os << indent << "DataStorage " << this << " is managing " << all->Size() << " objects. List of objects:" << std::endl; for (mitk::DataStorage::SetOfObjects::ConstIterator allIt = all->Begin(); allIt != all->End(); allIt++) { std::string name; allIt.Value()->GetName(name); std::string datatype; if (allIt.Value()->GetData() != NULL) datatype = allIt.Value()->GetData()->GetNameOfClass(); os << indent << " " << allIt.Value().GetPointer() << "<" << datatype << ">: " << name << std::endl; mitk::DataStorage::SetOfObjects::ConstPointer parents = this->GetSources(allIt.Value()); if (parents->Size() > 0) { os << indent << " Direct sources: "; for (mitk::DataStorage::SetOfObjects::ConstIterator parentIt = parents->Begin(); parentIt != parents->End(); parentIt++) os << parentIt.Value().GetPointer() << ", "; os << std::endl; } mitk::DataStorage::SetOfObjects::ConstPointer derivations = this->GetDerivations(allIt.Value()); if (derivations->Size() > 0) { os << indent << " Direct derivations: "; for (mitk::DataStorage::SetOfObjects::ConstIterator derivationIt = derivations->Begin(); derivationIt != derivations->End(); derivationIt++) os << derivationIt.Value().GetPointer() << ", "; os << std::endl; } } os << std::endl; } mitk::DataStorage::SetOfObjects::ConstPointer mitk::DataStorage::FilterSetOfObjects( const SetOfObjects *set, const NodePredicateBase *condition) const { if (set == NULL) return NULL; mitk::DataStorage::SetOfObjects::Pointer result = mitk::DataStorage::SetOfObjects::New(); for (mitk::DataStorage::SetOfObjects::ConstIterator it = set->Begin(); it != set->End(); it++) if (condition == NULL || condition->CheckNode(it.Value()) == true) // alway copy the set, otherwise the iterator in mitk::DataStorage::Remove() will crash result->InsertElement(result->Size(), it.Value()); return mitk::DataStorage::SetOfObjects::ConstPointer(result); } const mitk::DataNode::GroupTagList mitk::DataStorage::GetGroupTags() const { DataNode::GroupTagList result; SetOfObjects::ConstPointer all = this->GetAll(); if (all.IsNull()) return result; for (mitk::DataStorage::SetOfObjects::ConstIterator nodeIt = all->Begin(); nodeIt != all->End(); nodeIt++) // for each node { mitk::PropertyList *pl = nodeIt.Value()->GetPropertyList(); for (mitk::PropertyList::PropertyMap::const_iterator propIt = pl->GetMap()->begin(); propIt != pl->GetMap()->end(); ++propIt) if (dynamic_cast(propIt->second.GetPointer()) != NULL) result.insert(propIt->first); } return result; } void mitk::DataStorage::EmitAddNodeEvent(const mitk::DataNode *node) { AddNodeEvent.Send(node); } void mitk::DataStorage::EmitRemoveNodeEvent(const mitk::DataNode *node) { RemoveNodeEvent.Send(node); } void mitk::DataStorage::OnNodeInteractorChanged(itk::Object *caller, const itk::EventObject &) { const mitk::DataNode *_Node = dynamic_cast(caller); if (_Node) { InteractorChangedNodeEvent.Send(_Node); } } void mitk::DataStorage::OnNodeModifiedOrDeleted(const itk::Object *caller, const itk::EventObject &event) { if (m_BlockNodeModifiedEvents) return; const mitk::DataNode *_Node = dynamic_cast(caller); if (_Node) { const itk::ModifiedEvent *modEvent = dynamic_cast(&event); if (modEvent) ChangedNodeEvent.Send(_Node); else DeleteNodeEvent.Send(_Node); } } void mitk::DataStorage::AddListeners(const mitk::DataNode *_Node) { itk::MutexLockHolder locked(m_MutexOne); // node must not be 0 and must not be yet registered mitk::DataNode *NonConstNode = const_cast(_Node); if (_Node && m_NodeModifiedObserverTags.find(NonConstNode) == m_NodeModifiedObserverTags.end()) { itk::MemberCommand::Pointer nodeModifiedCommand = itk::MemberCommand::New(); nodeModifiedCommand->SetCallbackFunction(this, &mitk::DataStorage::OnNodeModifiedOrDeleted); m_NodeModifiedObserverTags[NonConstNode] = NonConstNode->AddObserver(itk::ModifiedEvent(), nodeModifiedCommand); itk::MemberCommand::Pointer interactorChangedCommand = itk::MemberCommand::New(); interactorChangedCommand->SetCallbackFunction(this, &mitk::DataStorage::OnNodeInteractorChanged); m_NodeInteractorChangedObserverTags[NonConstNode] = NonConstNode->AddObserver(mitk::DataNode::InteractorChangedEvent(), interactorChangedCommand); // add itk delete listener on datastorage itk::MemberCommand::Pointer deleteCommand = itk::MemberCommand::New(); deleteCommand->SetCallbackFunction(this, &mitk::DataStorage::OnNodeModifiedOrDeleted); // add observer m_NodeDeleteObserverTags[NonConstNode] = NonConstNode->AddObserver(itk::DeleteEvent(), deleteCommand); } } void mitk::DataStorage::RemoveListeners(const mitk::DataNode *_Node) { itk::MutexLockHolder locked(m_MutexOne); // node must not be 0 and must be registered mitk::DataNode *NonConstNode = const_cast(_Node); if (_Node && m_NodeModifiedObserverTags.find(NonConstNode) != m_NodeModifiedObserverTags.end()) { // const cast is bad! but sometimes it is necessary. removing an observer does not really // touch the internal state NonConstNode->RemoveObserver(m_NodeModifiedObserverTags.find(NonConstNode)->second); NonConstNode->RemoveObserver(m_NodeDeleteObserverTags.find(NonConstNode)->second); NonConstNode->RemoveObserver(m_NodeInteractorChangedObserverTags.find(NonConstNode)->second); m_NodeModifiedObserverTags.erase(NonConstNode); m_NodeDeleteObserverTags.erase(NonConstNode); m_NodeInteractorChangedObserverTags.erase(NonConstNode); } } mitk::TimeGeometry::Pointer mitk::DataStorage::ComputeBoundingGeometry3D(const SetOfObjects *input, const char *boolPropertyKey, const mitk::BaseRenderer *renderer, const char *boolPropertyKey2) const { if (input == NULL) throw std::invalid_argument("DataStorage: input is invalid"); BoundingBox::PointsContainer::Pointer pointscontainer = BoundingBox::PointsContainer::New(); BoundingBox::PointIdentifier pointid = 0; Point3D point; Vector3D minSpacing; minSpacing.Fill(itk::NumericTraits::max()); ScalarType stmin, stmax; stmin = itk::NumericTraits::NonpositiveMin(); stmax = itk::NumericTraits::max(); ScalarType minimalIntervallSize = stmax; ScalarType minimalTime = stmax; ScalarType maximalTime = 0; // Needed for check of zero bounding boxes mitk::ScalarType nullpoint[] = {0, 0, 0, 0, 0, 0}; BoundingBox::BoundsArrayType itkBoundsZero(nullpoint); for (SetOfObjects::ConstIterator it = input->Begin(); it != input->End(); ++it) { DataNode::Pointer node = it->Value(); if ((node.IsNotNull()) && (node->GetData() != NULL) && (node->GetData()->IsEmpty() == false) && node->IsOn(boolPropertyKey, renderer) && node->IsOn(boolPropertyKey2, renderer)) { const TimeGeometry *timeGeometry = node->GetData()->GetUpdatedTimeGeometry(); if (timeGeometry != NULL) { // bounding box (only if non-zero) BoundingBox::BoundsArrayType itkBounds = timeGeometry->GetBoundingBoxInWorld()->GetBounds(); if (itkBounds == itkBoundsZero) { continue; } unsigned char i; for (i = 0; i < 8; ++i) { point = timeGeometry->GetCornerPointInWorld(i); if (point[0] * point[0] + point[1] * point[1] + point[2] * point[2] < large) pointscontainer->InsertElement(pointid++, point); else { itkGenericOutputMacro(<< "Unrealistically distant corner point encountered. Ignored. Node: " << node); } } try { // time bounds // iterate over all time steps // Attention: Objects with zero bounding box are not respected in time bound calculation for (TimeStepType i = 0; i < timeGeometry->CountTimeSteps(); i++) { - Vector3D spacing = node->GetData()->GetGeometry(i)->GetSpacing(); - for (int axis = 0; axis < 3; ++axis) + // We must not use 'node->GetData()->GetGeometry(i)->GetSpacing()' here, as it returns the spacing + // in its original space, which, in case of an image geometry, can have the values in different + // order than in world space. For the further calculations, we need to have the spacing values + // in world coordinate order (sag-cor-ax). + Vector3D spacing; + spacing.Fill(1.0); + node->GetData()->GetGeometry(i)->IndexToWorld(spacing, spacing); + for (int axis = 0; axis < 3; ++ axis) { - if (spacing[axis] < minSpacing[axis]) - minSpacing[axis] = spacing[axis]; + ScalarType space = std::abs(spacing[axis]); + if (space < minSpacing[axis]) + { + minSpacing[axis] = space; + } } const TimeBounds &curTimeBounds = node->GetData()->GetTimeGeometry()->GetTimeBounds(i); // get the minimal time of all objects in the DataStorage if ((curTimeBounds[0] < minimalTime) && (curTimeBounds[0] > stmin)) { minimalTime = curTimeBounds[0]; } // get the maximal time of all objects in the DataStorage if ((curTimeBounds[1] > maximalTime) && (curTimeBounds[1] < stmax)) { maximalTime = curTimeBounds[1]; } // get the minimal TimeBound of all time steps of the current DataNode if (curTimeBounds[1] - curTimeBounds[0] < minimalIntervallSize) { minimalIntervallSize = curTimeBounds[1] - curTimeBounds[0]; } } } catch (itk::ExceptionObject &e) { MITK_ERROR << e << std::endl; } } } } BoundingBox::Pointer result = BoundingBox::New(); result->SetPoints(pointscontainer); result->ComputeBoundingBox(); // compute the number of time steps unsigned int numberOfTimeSteps = 1; if (maximalTime == 0) // make sure that there is at least one time sliced geometry in the data storage { minimalTime = 0; maximalTime = 1; minimalIntervallSize = 1; } numberOfTimeSteps = static_cast((maximalTime - minimalTime) / minimalIntervallSize); TimeGeometry::Pointer timeGeometry = NULL; if (result->GetPoints()->Size() > 0) { // Initialize a geometry of a single time step Geometry3D::Pointer geometry = Geometry3D::New(); geometry->Initialize(); // correct bounding-box (is now in mm, should be in index-coordinates) // according to spacing BoundingBox::BoundsArrayType bounds = result->GetBounds(); - int i; - for (i = 0; i < 6; ++i) + AffineTransform3D::OutputVectorType offset; + for (int i = 0; i < 3; ++i) { - bounds[i] /= minSpacing[i / 2]; + offset[i] = bounds[i * 2]; + bounds[i * 2] = 0.0; + bounds[i * 2 + 1] = (bounds[i * 2 + 1] - offset[i]) / minSpacing[i]; } + geometry->GetIndexToWorldTransform()->SetOffset(offset); geometry->SetBounds(bounds); geometry->SetSpacing(minSpacing); // Initialize the time sliced geometry timeGeometry = ProportionalTimeGeometry::New(); dynamic_cast(timeGeometry.GetPointer())->Initialize(geometry, numberOfTimeSteps); dynamic_cast(timeGeometry.GetPointer())->SetFirstTimePoint(minimalTime); dynamic_cast(timeGeometry.GetPointer())->SetStepDuration(minimalIntervallSize); } return timeGeometry; } mitk::TimeGeometry::Pointer mitk::DataStorage::ComputeBoundingGeometry3D(const char *boolPropertyKey, const mitk::BaseRenderer *renderer, const char *boolPropertyKey2) const { return this->ComputeBoundingGeometry3D(this->GetAll(), boolPropertyKey, renderer, boolPropertyKey2); } mitk::TimeGeometry::Pointer mitk::DataStorage::ComputeVisibleBoundingGeometry3D(const mitk::BaseRenderer *renderer, const char *boolPropertyKey) { return ComputeBoundingGeometry3D("visible", renderer, boolPropertyKey); } mitk::BoundingBox::Pointer mitk::DataStorage::ComputeBoundingBox(const char *boolPropertyKey, const mitk::BaseRenderer *renderer, const char *boolPropertyKey2) { BoundingBox::PointsContainer::Pointer pointscontainer = BoundingBox::PointsContainer::New(); BoundingBox::PointIdentifier pointid = 0; Point3D point; // Needed for check of zero bounding boxes mitk::ScalarType nullpoint[] = {0, 0, 0, 0, 0, 0}; BoundingBox::BoundsArrayType itkBoundsZero(nullpoint); SetOfObjects::ConstPointer all = this->GetAll(); for (SetOfObjects::ConstIterator it = all->Begin(); it != all->End(); ++it) { DataNode::Pointer node = it->Value(); if ((node.IsNotNull()) && (node->GetData() != NULL) && (node->GetData()->IsEmpty() == false) && node->IsOn(boolPropertyKey, renderer) && node->IsOn(boolPropertyKey2, renderer)) { const TimeGeometry *geometry = node->GetData()->GetUpdatedTimeGeometry(); if (geometry != NULL) { // bounding box (only if non-zero) BoundingBox::BoundsArrayType itkBounds = geometry->GetBoundingBoxInWorld()->GetBounds(); if (itkBounds == itkBoundsZero) { continue; } unsigned char i; for (i = 0; i < 8; ++i) { point = geometry->GetCornerPointInWorld(i); if (point[0] * point[0] + point[1] * point[1] + point[2] * point[2] < large) pointscontainer->InsertElement(pointid++, point); else { itkGenericOutputMacro(<< "Unrealistically distant corner point encountered. Ignored. Node: " << node); } } } } } BoundingBox::Pointer result = BoundingBox::New(); result->SetPoints(pointscontainer); result->ComputeBoundingBox(); return result; } mitk::TimeBounds mitk::DataStorage::ComputeTimeBounds(const char *boolPropertyKey, const mitk::BaseRenderer *renderer, const char *boolPropertyKey2) { TimeBounds timeBounds; ScalarType stmin, stmax, cur; stmin = itk::NumericTraits::NonpositiveMin(); stmax = itk::NumericTraits::max(); timeBounds[0] = stmax; timeBounds[1] = stmin; SetOfObjects::ConstPointer all = this->GetAll(); for (SetOfObjects::ConstIterator it = all->Begin(); it != all->End(); ++it) { DataNode::Pointer node = it->Value(); if ((node.IsNotNull()) && (node->GetData() != NULL) && (node->GetData()->IsEmpty() == false) && node->IsOn(boolPropertyKey, renderer) && node->IsOn(boolPropertyKey2, renderer)) { const TimeGeometry *geometry = node->GetData()->GetUpdatedTimeGeometry(); if (geometry != NULL) { const TimeBounds &curTimeBounds = geometry->GetTimeBounds(); cur = curTimeBounds[0]; // is it after -infinity, but before everything else that we found until now? if ((cur > stmin) && (cur < timeBounds[0])) timeBounds[0] = cur; cur = curTimeBounds[1]; // is it before infinity, but after everything else that we found until now? if ((cur < stmax) && (cur > timeBounds[1])) timeBounds[1] = cur; } } } if (!(timeBounds[0] < stmax)) { timeBounds[0] = stmin; timeBounds[1] = stmax; } return timeBounds; } void mitk::DataStorage::BlockNodeModifiedEvents(bool block) { m_BlockNodeModifiedEvents = block; } diff --git a/Modules/DiffusionImaging/DiffusionCore/cmdapps/CMakeLists.txt b/Modules/DiffusionImaging/DiffusionCore/cmdapps/CMakeLists.txt index 11f0c2aa2b..7fbf48d40b 100644 --- a/Modules/DiffusionImaging/DiffusionCore/cmdapps/CMakeLists.txt +++ b/Modules/DiffusionImaging/DiffusionCore/cmdapps/CMakeLists.txt @@ -1,44 +1,44 @@ option(BUILD_DiffusionMiniApps "Build commandline tools for diffusion" OFF) if(BUILD_DiffusionMiniApps OR MITK_BUILD_ALL_APPS) # needed include directories include_directories( ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ) # list of diffusion miniapps # if an app requires additional dependencies # they are added after a "^^" and separated by "_" set( diffusionminiapps DwiDenoising^^ ImageResampler^^ ExportShImage^^ CopyGeometry^^ DiffusionIndices^^ QballReconstruction^^ Registration^^ TensorReconstruction^^ TensorDerivedMapsExtraction^^ DiffusionDICOMLoader^^ - DiffusionIVIMFit^^ + DiffusionKurtosisFit^^ ) foreach(diffusionminiapp ${diffusionminiapps}) # extract mini app name and dependencies string(REPLACE "^^" "\\;" miniapp_info ${diffusionminiapp}) set(miniapp_info_list ${miniapp_info}) list(GET miniapp_info_list 0 appname) list(GET miniapp_info_list 1 raw_dependencies) string(REPLACE "_" "\\;" dependencies "${raw_dependencies}") set(dependencies_list ${dependencies}) mitkFunctionCreateCommandLineApp( NAME ${appname} DEPENDS MitkCore MitkDiffusionCore ${dependencies_list} PACKAGE_DEPENDS ITK ) endforeach() endif() diff --git a/Modules/DiffusionImaging/DiffusionCore/cmdapps/DiffusionIVIMFit.cpp b/Modules/DiffusionImaging/DiffusionCore/cmdapps/DiffusionKurtosisFit.cpp similarity index 60% rename from Modules/DiffusionImaging/DiffusionCore/cmdapps/DiffusionIVIMFit.cpp rename to Modules/DiffusionImaging/DiffusionCore/cmdapps/DiffusionKurtosisFit.cpp index 6c49d7a7d0..d1886d0469 100644 --- a/Modules/DiffusionImaging/DiffusionCore/cmdapps/DiffusionIVIMFit.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/cmdapps/DiffusionKurtosisFit.cpp @@ -1,207 +1,253 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include #include +#include "mitkImage.h" +#include +#include +#include +#include "mitkIOUtil.h" + +#include +#include + + //vnl_includes #include "vnl/vnl_math.h" #include "vnl/vnl_cost_function.h" #include "vnl/vnl_least_squares_function.h" #include "vnl/algo/vnl_lbfgsb.h" #include "vnl/algo/vnl_lbfgs.h" #include "vnl/algo/vnl_levenberg_marquardt.h" typedef mitk::DiffusionPropertyHelper DPH; #include #include #include #include #include #include + + DPH::ImageType::Pointer GetBlurredVectorImage( DPH::ImageType::Pointer vectorImage, double sigma) { typedef itk::DiscreteGaussianImageFilter< itk::Image, itk::Image > GaussianFilterType; typedef itk::VectorIndexSelectionCastImageFilter< DPH::ImageType, itk::Image > IndexSelectionType; IndexSelectionType::Pointer indexSelectionFilter = IndexSelectionType::New(); indexSelectionFilter->SetInput( vectorImage ); typedef itk::ComposeImageFilter< itk::Image, DPH::ImageType > ComposeFilterType; ComposeFilterType::Pointer vec_composer = ComposeFilterType::New(); for( unsigned int i=0; iGetVectorLength(); ++i) { GaussianFilterType::Pointer gaussian_filter = GaussianFilterType::New(); indexSelectionFilter->SetIndex( i ); gaussian_filter->SetInput( indexSelectionFilter->GetOutput() ); gaussian_filter->SetVariance( sigma ); vec_composer->SetInput(i, gaussian_filter->GetOutput() ); gaussian_filter->Update(); } try { vec_composer->Update(); } catch(const itk::ExceptionObject &e) { mitkThrow() << "[VectorImage.GaussianSmoothing] !! Failed with ITK Exception: " << e.what(); } DPH::ImageType::Pointer smoothed_vector = vec_composer->GetOutput(); -/* + + /* itk::ImageFileWriter< DPH::ImageType >::Pointer writer = itk::ImageFileWriter< DPH::ImageType >::New(); writer->SetInput( smoothed_vector ); writer->SetFileName( "/tmp/itk_smoothed_vector.nrrd"); writer->Update();*/ return smoothed_vector; } -void KurtosisMapComputation( mitk::Image::Pointer input, std::string output_prefix ) +void KurtosisMapComputation( mitk::Image::Pointer input, + std::string output_prefix , + std::string output_type, + std::string maskPath, + bool omitBZero, + double lower, + double upper ) { DPH::ImageType::Pointer vectorImage = DPH::ImageType::New(); mitk::CastToItkImage( input, vectorImage ); - - typedef itk::DiffusionKurtosisReconstructionImageFilter< short, double > KurtosisFilterType; KurtosisFilterType::Pointer kurtosis_filter = KurtosisFilterType::New(); - kurtosis_filter->SetInput( GetBlurredVectorImage( vectorImage, 1.5 ) ); kurtosis_filter->SetReferenceBValue( DPH::GetReferenceBValue( input.GetPointer() ) ); kurtosis_filter->SetGradientDirections( DPH::GetGradientContainer( input.GetPointer() ) ); - kurtosis_filter->SetNumberOfThreads(1); +// kurtosis_filter->SetNumberOfThreads(1); + kurtosis_filter->SetOmitUnweightedValue(omitBZero); + kurtosis_filter->SetBoundariesForKurtosis(-lower,upper); +// kurtosis_filter->SetInitialSolution(const vnl_vector& x0 ); - KurtosisFilterType::OutputImageRegionType o_region; - KurtosisFilterType::OutputImageRegionType::SizeType o_size; - KurtosisFilterType::OutputImageRegionType::IndexType o_index; - o_index.Fill(0); o_size.Fill(0); - o_index[0] = 48; o_index[1] = 18; o_index[2] = 12; - o_size[0] = 16; o_size[1] = 24; o_size[2] = 1; - - o_region.SetSize( o_size ); - o_region.SetIndex( o_index ); - kurtosis_filter->SetMapOutputRegion( o_region ); + if(maskPath != "") + { + mitk::Image::Pointer segmentation; + segmentation = mitk::IOUtil::LoadImage(maskPath); + typedef itk::Image< short , 3> MaskImageType; + MaskImageType::Pointer vectorSeg = MaskImageType::New() ; + mitk::CastToItkImage( segmentation, vectorSeg ); + kurtosis_filter->SetImageMask(vectorSeg) ; + } try { kurtosis_filter->Update(); } catch( const itk::ExceptionObject& e) { mitkThrow() << "Kurtosis fit failed with an ITK Exception: " << e.what(); } mitk::Image::Pointer d_image = mitk::Image::New(); d_image->InitializeByItk( kurtosis_filter->GetOutput(0) ); d_image->SetVolume( kurtosis_filter->GetOutput(0)->GetBufferPointer() ); mitk::Image::Pointer k_image = mitk::Image::New(); k_image->InitializeByItk( kurtosis_filter->GetOutput(1) ); k_image->SetVolume( kurtosis_filter->GetOutput(1)->GetBufferPointer() ); - std::string outputD_FileName = output_prefix + "_ADC_map.nrrd"; - std::string outputK_FileName = output_prefix + "_AKC_map.nrrd"; + std::string outputD_FileName = output_prefix + "_ADC_map." + output_type; + std::string outputK_FileName = output_prefix + "_AKC_map." + output_type; try { - mitk::IOUtil::Save( d_image, outputD_FileName ); - mitk::IOUtil::Save( k_image, outputK_FileName ); + mitk::IOUtil::Save( d_image, outputD_FileName ); + mitk::IOUtil::Save( k_image, outputK_FileName ); } catch( const itk::ExceptionObject& e) { mitkThrow() << "Failed to save the KurtosisFit Results due to exception: " << e.what(); } } int main( int argc, char* argv[] ) { - mitkCommandLineParser parser; - parser.setTitle("Diffusion IVIM (Kurtosis) Fit"); + parser.setTitle("Diffusion Kurtosis Fit"); parser.setCategory("Signal Reconstruction"); parser.setContributor("MIC"); - parser.setDescription("Fitting of IVIM / Kurtosis"); - + parser.setDescription("Fitting Kurtosis"); parser.setArgumentPrefix("--","-"); + // mandatory arguments parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input: ", "input image (DWI)", us::Any(), false); parser.addArgument("output", "o", mitkCommandLineParser::String, "Output Preifx: ", "Prefix for the output images, will append _f, _K, _D accordingly ", us::Any(), false); - parser.addArgument("fit", "f", mitkCommandLineParser::String, "Input: ", "input image (DWI)", us::Any(), false); + parser.addArgument("output_type", "ot", mitkCommandLineParser::String, "Output Type: ", "choose data type of output image, e.g. '.nii' or '.nrrd' ", us::Any(), false); // optional arguments - parser.addArgument("mask", "m", mitkCommandLineParser::InputFile, "Masking Image: ", "ROI (segmentation)", us::Any(), true); + parser.addArgument("mask", "m", mitkCommandLineParser::InputFile, "Masking Image: ", "ROI (segmentation)", us::Any()); + parser.addArgument("help", "h", mitkCommandLineParser::Bool, "Help", "Show this help text"); + parser.addArgument("omitbzero", "om", mitkCommandLineParser::Bool, "Omit b0:", "Omit b0 value during fit (default = false)", us::Any()); + parser.addArgument("lowerkbound", "kl", mitkCommandLineParser::Float, "lower Kbound:", "Set (unsigned) lower boundary for Kurtosis parameter (default = -1000)", us::Any()); + parser.addArgument("upperkbound", "ku", mitkCommandLineParser::Float, "upper Kbound:", "Set upper boundary for Kurtosis parameter (default = 1000)", us::Any()); std::map parsedArgs = parser.parseArguments(argc, argv); - if (parsedArgs.size()==0) - return EXIT_FAILURE; + + if (parsedArgs.size()==0 || parsedArgs.count("help") || parsedArgs.count("h")){ + std::cout << parser.helpText(); + return EXIT_SUCCESS; + } // mandatory arguments std::string inFileName = us::any_cast(parsedArgs["input"]); std::string out_prefix = us::any_cast(parsedArgs["output"]); - std::string fit_name = us::any_cast(parsedArgs["fit"]); + std::string maskPath = ""; mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inFileName); - if( !DPH::IsDiffusionWeightedImage( inputImage ) ) + + bool omitBZero = false; + double lower = -1000; + double upper = 1000; + std::string out_type = "nrrd"; + + if (parsedArgs.count("mask") || parsedArgs.count("m")) { - MITK_ERROR("DiffusionIVIMFit.Input") << "No valid diffusion-weighted image provided, failed to load " << inFileName << " as DW Image. Aborting..."; - return EXIT_FAILURE; + maskPath = us::any_cast(parsedArgs["mask"]); } - if( fit_name == "Kurtosis" ) + if (parsedArgs.count("output_type") || parsedArgs.count("ot")) { - MITK_INFO("DiffusionIVIMFit.Main") << "-----[ Kurtosis Fit ]-----"; + out_type = us::any_cast(parsedArgs["output_type"]); + } - KurtosisMapComputation( inputImage, out_prefix ); + if (parsedArgs.count("omitbzero") || parsedArgs.count("om")) + { + omitBZero = us::any_cast(parsedArgs["omitbzero"]); + } + if (parsedArgs.count("lowerkbound") || parsedArgs.count("kl")) + { + lower = us::any_cast(parsedArgs["lowerkbound"]); } - else if (fit_name == "IVIM" ) + + if (parsedArgs.count("upperkbound") || parsedArgs.count("ku")) { - MITK_INFO("DiffusionIVIMFit.Main") << "IVIM Fit not fully implemented yet. Aborting..."; - return EXIT_FAILURE; + upper = us::any_cast(parsedArgs["upperkbound"]); } - else + + if( !DPH::IsDiffusionWeightedImage( inputImage ) ) { - MITK_ERROR("DiffusionIVIMFit.Main") << "Unrecognized option: " << fit_name << ". Valid values [\"IVIM\", \"Kurtosis\"] \n Aborting... \n"; + MITK_ERROR("DiffusionIVIMFit.Input") << "No valid diffusion-weighted image provided, failed to load " << inFileName << " as DW Image. Aborting..."; return EXIT_FAILURE; - } + + +KurtosisMapComputation( inputImage, + out_prefix , + out_type, + maskPath, + omitBZero, + lower, + upper); + }