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<mitk::GroupTagProperty *>(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<const mitk::DataNode *>(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<const mitk::DataNode *>(caller);
if (_Node)
{
const itk::ModifiedEvent *modEvent = dynamic_cast<const itk::ModifiedEvent *>(&event);
if (modEvent)
ChangedNodeEvent.Send(_Node);
else
DeleteNodeEvent.Send(_Node);
}
}
void mitk::DataStorage::AddListeners(const mitk::DataNode *_Node)
{
itk::MutexLockHolder<itk::SimpleFastMutexLock> locked(m_MutexOne);
// node must not be 0 and must not be yet registered
mitk::DataNode *NonConstNode = const_cast<mitk::DataNode *>(_Node);
if (_Node && m_NodeModifiedObserverTags.find(NonConstNode) == m_NodeModifiedObserverTags.end())
{
itk::MemberCommand<mitk::DataStorage>::Pointer nodeModifiedCommand = itk::MemberCommand<mitk::DataStorage>::New();
nodeModifiedCommand->SetCallbackFunction(this, &mitk::DataStorage::OnNodeModifiedOrDeleted);
m_NodeModifiedObserverTags[NonConstNode] = NonConstNode->AddObserver(itk::ModifiedEvent(), nodeModifiedCommand);
itk::MemberCommand<mitk::DataStorage>::Pointer interactorChangedCommand =
itk::MemberCommand<mitk::DataStorage>::New();
interactorChangedCommand->SetCallbackFunction(this, &mitk::DataStorage::OnNodeInteractorChanged);
m_NodeInteractorChangedObserverTags[NonConstNode] =
NonConstNode->AddObserver(mitk::DataNode::InteractorChangedEvent(), interactorChangedCommand);
// add itk delete listener on datastorage
itk::MemberCommand<mitk::DataStorage>::Pointer deleteCommand = itk::MemberCommand<mitk::DataStorage>::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<itk::SimpleFastMutexLock> locked(m_MutexOne);
// node must not be 0 and must be registered
mitk::DataNode *NonConstNode = const_cast<mitk::DataNode *>(_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<mitk::ScalarType>::max());
ScalarType stmin, stmax;
stmin = itk::NumericTraits<mitk::ScalarType>::NonpositiveMin();
stmax = itk::NumericTraits<mitk::ScalarType>::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<unsigned int>((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<ProportionalTimeGeometry *>(timeGeometry.GetPointer())->Initialize(geometry, numberOfTimeSteps);
dynamic_cast<ProportionalTimeGeometry *>(timeGeometry.GetPointer())->SetFirstTimePoint(minimalTime);
dynamic_cast<ProportionalTimeGeometry *>(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<mitk::ScalarType>::NonpositiveMin();
stmax = itk::NumericTraits<mitk::ScalarType>::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 <mitkIOUtil.h>
#include <mitkDiffusionPropertyHelper.h>
#include <mitkImageCaster.h>
#include <itkDiffusionKurtosisReconstructionImageFilter.h>
+#include "mitkImage.h"
+#include <iostream>
+#include <usAny.h>
+#include <fstream>
+#include "mitkIOUtil.h"
+
+#include <itkFileTools.h>
+#include <itksys/SystemTools.hxx>
+
+
//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 <itkGaussianBlurImageFunction.h>
#include <itkUnaryFunctorImageFilter.h>
#include <itkImageFileWriter.h>
#include <itkVectorIndexSelectionCastImageFilter.h>
#include <itkComposeImageFilter.h>
#include <itkDiscreteGaussianImageFilter.h>
+
+
DPH::ImageType::Pointer GetBlurredVectorImage( DPH::ImageType::Pointer vectorImage, double sigma)
{
typedef itk::DiscreteGaussianImageFilter< itk::Image<DPH::DiffusionPixelType, 3 >, itk::Image<DPH::DiffusionPixelType, 3 > > GaussianFilterType;
typedef itk::VectorIndexSelectionCastImageFilter< DPH::ImageType, itk::Image<DPH::DiffusionPixelType, 3 > > IndexSelectionType;
IndexSelectionType::Pointer indexSelectionFilter = IndexSelectionType::New();
indexSelectionFilter->SetInput( vectorImage );
typedef itk::ComposeImageFilter< itk::Image<DPH::DiffusionPixelType, 3>, DPH::ImageType > ComposeFilterType;
ComposeFilterType::Pointer vec_composer = ComposeFilterType::New();
for( unsigned int i=0; i<vectorImage->GetVectorLength(); ++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<double>& 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<std::string, us::Any> 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<std::string>(parsedArgs["input"]);
std::string out_prefix = us::any_cast<std::string>(parsedArgs["output"]);
- std::string fit_name = us::any_cast<std::string>(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<std::string>(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<std::string>(parsedArgs["output_type"]);
+ }
- KurtosisMapComputation( inputImage, out_prefix );
+ if (parsedArgs.count("omitbzero") || parsedArgs.count("om"))
+ {
+ omitBZero = us::any_cast<bool>(parsedArgs["omitbzero"]);
+ }
+ if (parsedArgs.count("lowerkbound") || parsedArgs.count("kl"))
+ {
+ lower = us::any_cast<float>(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<float>(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);
+
}