diff --git a/Modules/CEST/files.cmake b/Modules/CEST/files.cmake
index 34a71a20af..e140e702c8 100644
--- a/Modules/CEST/files.cmake
+++ b/Modules/CEST/files.cmake
@@ -1,13 +1,14 @@
file(GLOB_RECURSE H_FILES RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" "${CMAKE_CURRENT_SOURCE_DIR}/include/*")
set(CPP_FILES
mitkCESTImageNormalizationFilter.cpp
mitkCustomTagParser.cpp
mitkCESTImageDetectionHelper.cpp
+ mitkExtractCESTOffset.cpp
)
set(RESOURCE_FILES
1416.json
1485.json
1494.json
)
diff --git a/Modules/CEST/include/mitkExtractCESTOffset.h b/Modules/CEST/include/mitkExtractCESTOffset.h
new file mode 100644
index 0000000000..6434060a8b
--- /dev/null
+++ b/Modules/CEST/include/mitkExtractCESTOffset.h
@@ -0,0 +1,39 @@
+/*============================================================================
+
+The Medical Imaging Interaction Toolkit (MITK)
+
+Copyright (c) German Cancer Research Center (DKFZ)
+All rights reserved.
+
+Use of this source code is governed by a 3-clause BSD license that can be
+found in the LICENSE file.
+
+============================================================================*/
+
+#ifndef __MITK_EXTRACT_CEST_OFFSET_H_
+#define __MITK_EXTRACT_CEST_OFFSET_H_
+
+#include <mitkBaseData.h>
+
+#include "MitkCESTExports.h"
+
+namespace mitk
+{
+ /**Helper function that gets the CEST offset property ("CEST.Offsets") from the input
+ image as vector of ScalarType.
+ If it is not possible to generate/get the offset an mitk::Excpetion will be thrown.
+ The values of the vector are in [ppm].
+ @post Number of extracted offsets equal the number of timesteps of the image.
+ */
+ MITKCEST_EXPORT std::vector<ScalarType> ExtractCESTOffset(const BaseData* image);
+
+ /**Helper function that gets the CEST offset property ("CEST.TREC") from the input image as vector of ScalarType.
+ If it is not possible to generate/get the T1 times an mitk::Excpetion will be thrown.
+ The values of the vector are in [sec]. In the property they are stored in [ms] and scaled appropriately
+ before returning.
+ @post Number of extracted T1 times equal the number of timesteps of the image.
+ */
+ MITKCEST_EXPORT std::vector<ScalarType> ExtractCESTT1Time(const BaseData* image);
+}
+
+#endif
diff --git a/Modules/CEST/src/mitkCESTImageNormalizationFilter.cpp b/Modules/CEST/src/mitkCESTImageNormalizationFilter.cpp
index 8648ae9bc3..b6f1b1f1ad 100644
--- a/Modules/CEST/src/mitkCESTImageNormalizationFilter.cpp
+++ b/Modules/CEST/src/mitkCESTImageNormalizationFilter.cpp
@@ -1,229 +1,202 @@
/*============================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center (DKFZ)
All rights reserved.
Use of this source code is governed by a 3-clause BSD license that can be
found in the LICENSE file.
============================================================================*/
#include "mitkCESTImageNormalizationFilter.h"
#include <mitkCustomTagParser.h>
+#include <mitkExtractCESTOffset.h>
#include <mitkImage.h>
#include <mitkImageAccessByItk.h>
#include <mitkImageCast.h>
#include <boost/algorithm/string.hpp>
mitk::CESTImageNormalizationFilter::CESTImageNormalizationFilter()
{
}
mitk::CESTImageNormalizationFilter::~CESTImageNormalizationFilter()
{
}
void mitk::CESTImageNormalizationFilter::GenerateData()
{
mitk::Image::ConstPointer inputImage = this->GetInput(0);
if ((inputImage->GetDimension() != 4))
{
mitkThrow() << "mitk::CESTImageNormalizationFilter:GenerateData works only with 4D images, sorry.";
return;
}
auto resultMitkImage = this->GetOutput();
AccessFixedDimensionByItk(inputImage, NormalizeTimeSteps, 4);
auto originalTimeGeometry = this->GetInput()->GetTimeGeometry();
auto resultTimeGeometry = mitk::ProportionalTimeGeometry::New();
unsigned int numberOfNonM0s = m_NonM0Indices.size();
resultTimeGeometry->Expand(numberOfNonM0s);
for (unsigned int index = 0; index < numberOfNonM0s; ++index)
{
resultTimeGeometry->SetTimeStepGeometry(originalTimeGeometry->GetGeometryCloneForTimeStep(m_NonM0Indices.at(index)), index);
}
resultMitkImage->SetTimeGeometry(resultTimeGeometry);
resultMitkImage->SetPropertyList(this->GetInput()->GetPropertyList()->Clone());
resultMitkImage->GetPropertyList()->SetStringProperty(mitk::CEST_PROPERTY_NAME_OFFSETS().c_str(), m_RealOffsets.c_str());
// remove uids
resultMitkImage->GetPropertyList()->DeleteProperty("DICOM.0008.0018");
resultMitkImage->GetPropertyList()->DeleteProperty("DICOM.0020.000D");
resultMitkImage->GetPropertyList()->DeleteProperty("DICOM.0020.000E");
}
-std::vector<double> ExtractOffsets(const mitk::Image* image)
-{
- std::vector<double> result;
-
- if (image)
- {
- std::string offsets = "";
- std::vector<std::string> parts;
- if (image->GetPropertyList()->GetStringProperty(mitk::CEST_PROPERTY_NAME_OFFSETS().c_str(), offsets) && !offsets.empty())
- {
- boost::algorithm::trim(offsets);
- boost::split(parts, offsets, boost::is_any_of(" "));
-
- for (auto part : parts)
- {
- std::istringstream iss(part);
- iss.imbue(std::locale("C"));
- double d;
- iss >> d;
- result.push_back(d);
- }
- }
- }
-
- return result;
-}
-
-
template <typename TPixel, unsigned int VImageDimension>
void mitk::CESTImageNormalizationFilter::NormalizeTimeSteps(const itk::Image<TPixel, VImageDimension>* image)
{
typedef itk::Image<TPixel, VImageDimension> ImageType;
typedef itk::Image<double, VImageDimension> OutputImageType;
- auto offsets = ExtractOffsets(this->GetInput());
+ auto offsets = ExtractCESTOffset(this->GetInput());
// determine normalization images
std::vector<unsigned int> mZeroIndices;
std::stringstream offsetsWithoutM0;
offsetsWithoutM0.imbue(std::locale("C"));
m_NonM0Indices.clear();
for (unsigned int index = 0; index < offsets.size(); ++index)
{
if ((offsets.at(index) < -299) || (offsets.at(index) > 299))
{
mZeroIndices.push_back(index);
}
else
{
offsetsWithoutM0 << offsets.at(index) << " ";
m_NonM0Indices.push_back(index);
}
}
auto resultImage = OutputImageType::New();
typename ImageType::RegionType targetEntireRegion = image->GetLargestPossibleRegion();
targetEntireRegion.SetSize(3, m_NonM0Indices.size());
resultImage->SetRegions(targetEntireRegion);
resultImage->Allocate();
resultImage->FillBuffer(0);
unsigned int numberOfTimesteps = image->GetLargestPossibleRegion().GetSize(3);
typename ImageType::RegionType lowerMZeroRegion = image->GetLargestPossibleRegion();
lowerMZeroRegion.SetSize(3, 1);
typename ImageType::RegionType upperMZeroRegion = image->GetLargestPossibleRegion();
upperMZeroRegion.SetSize(3, 1);
typename ImageType::RegionType sourceRegion = image->GetLargestPossibleRegion();
sourceRegion.SetSize(3, 1);
typename OutputImageType::RegionType targetRegion = resultImage->GetLargestPossibleRegion();
targetRegion.SetSize(3, 1);
unsigned int targetTimestep = 0;
for (unsigned int sourceTimestep = 0; sourceTimestep < numberOfTimesteps; ++sourceTimestep)
{
unsigned int lowerMZeroIndex = mZeroIndices[0];
unsigned int upperMZeroIndex = mZeroIndices[0];
for (unsigned int loop = 0; loop < mZeroIndices.size(); ++loop)
{
if (mZeroIndices[loop] <= sourceTimestep)
{
lowerMZeroIndex = mZeroIndices[loop];
}
if (mZeroIndices[loop] > sourceTimestep)
{
upperMZeroIndex = mZeroIndices[loop];
break;
}
}
bool isMZero = (lowerMZeroIndex == sourceTimestep);
double weight = 0.0;
if (lowerMZeroIndex == upperMZeroIndex)
{
weight = 1.0;
}
else
{
weight = 1.0 - double(sourceTimestep - lowerMZeroIndex) / double(upperMZeroIndex - lowerMZeroIndex);
}
if (isMZero)
{
//do nothing
}
else
{
lowerMZeroRegion.SetIndex(3, lowerMZeroIndex);
upperMZeroRegion.SetIndex(3, upperMZeroIndex);
sourceRegion.SetIndex(3, sourceTimestep);
targetRegion.SetIndex(3, targetTimestep);
itk::ImageRegionConstIterator<ImageType> lowerMZeroIterator(image, lowerMZeroRegion);
itk::ImageRegionConstIterator<ImageType> upperMZeroIterator(image, upperMZeroRegion);
itk::ImageRegionConstIterator<ImageType> sourceIterator(image, sourceRegion);
itk::ImageRegionIterator<OutputImageType> targetIterator(resultImage.GetPointer(), targetRegion);
while (!sourceIterator.IsAtEnd())
{
double normalizationFactor = weight * lowerMZeroIterator.Get() + (1.0 - weight) * upperMZeroIterator.Get();
if (mitk::Equal(normalizationFactor, 0))
{
targetIterator.Set(0);
}
else
{
targetIterator.Set(double(sourceIterator.Get()) / normalizationFactor);
}
++lowerMZeroIterator;
++upperMZeroIterator;
++sourceIterator;
++targetIterator;
}
++targetTimestep;
}
}
// get Pointer to output image
mitk::Image::Pointer resultMitkImage = this->GetOutput();
// write into output image
mitk::CastToMitkImage<OutputImageType>(resultImage, resultMitkImage);
m_RealOffsets = offsetsWithoutM0.str();
}
void mitk::CESTImageNormalizationFilter::GenerateOutputInformation()
{
mitk::Image::ConstPointer input = this->GetInput();
mitk::Image::Pointer output = this->GetOutput();
itkDebugMacro(<< "GenerateOutputInformation()");
}
bool mitk::IsNotNormalizedCESTImage(const Image* cestImage)
{
- auto offsets = ExtractOffsets(cestImage);
+ auto offsets = ExtractCESTOffset(cestImage);
- for (auto offset : offsets)
+ for (const auto& offset : offsets)
{
if (offset < -299 || offset > 299)
{
return true;
}
}
return false;
};
diff --git a/Modules/CEST/src/mitkExtractCESTOffset.cpp b/Modules/CEST/src/mitkExtractCESTOffset.cpp
new file mode 100644
index 0000000000..f0b3878b0a
--- /dev/null
+++ b/Modules/CEST/src/mitkExtractCESTOffset.cpp
@@ -0,0 +1,86 @@
+/*===================================================================
+
+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 "mitkExtractCESTOffset.h"
+#include "mitkCustomTagParser.h"
+
+#include <iterator>
+#include <regex>
+
+
+std::vector<mitk::ScalarType> mitk::ExtractCESTT1Time(const BaseData* image)
+{
+ std::vector<ScalarType> result;
+
+ auto prop = image->GetProperty(CEST_PROPERTY_NAME_TREC().c_str());
+ if (prop.IsNotNull())
+ {
+ auto valueStr = prop->GetValueAsString();
+
+ std::istringstream iss;
+ iss.imbue(std::locale("C"));
+ iss.str(valueStr);
+ double d;
+
+ while (iss >> d)
+ {
+ if (!iss.fail())
+ {
+ result.emplace_back(d);
+ }
+ }
+
+ if (result.size() != image->GetTimeSteps()) mitkThrow() << "Cannot determine T1 times. Selected input has an property \"" << CEST_PROPERTY_NAME_TREC() << "\" that don't match the number of time steps of input.";
+
+ for (auto& value : result)
+ {
+ value *= 0.001;
+ }
+ }
+ else mitkThrow() << "Cannot determine T1 time grid (TREC). Selected input has no property \"" << CEST_PROPERTY_NAME_TREC() << "\"";
+
+ return result;
+}
+
+std::vector<mitk::ScalarType> mitk::ExtractCESTOffset(const BaseData* image)
+{
+ std::vector<ScalarType> result;
+
+ auto prop = image->GetProperty(CEST_PROPERTY_NAME_OFFSETS().c_str());
+ if (prop.IsNotNull())
+ {
+ auto valueStr = prop->GetValueAsString();
+
+ std::istringstream iss;
+ iss.imbue(std::locale("C"));
+ iss.str(valueStr);
+ double d;
+
+ while (iss >> d)
+ {
+ if (!iss.fail())
+ {
+ result.emplace_back(d);
+ }
+ }
+
+ if (result.size() != image->GetTimeSteps()) mitkThrow() << "Cannot determine offset. Selected input has an property \"" << CEST_PROPERTY_NAME_OFFSETS() << "\" that don't match the number of time steps of input.";
+
+ }
+ else mitkThrow() << "Cannot determine frequency. Selected input has no property \"" << CEST_PROPERTY_NAME_OFFSETS() << "\"";
+
+ return result;
+}
diff --git a/Plugins/org.mitk.gui.qt.cest/src/internal/QmitkCESTNormalizeView.cpp b/Plugins/org.mitk.gui.qt.cest/src/internal/QmitkCESTNormalizeView.cpp
index 05c8fb560a..d01c8204ca 100644
--- a/Plugins/org.mitk.gui.qt.cest/src/internal/QmitkCESTNormalizeView.cpp
+++ b/Plugins/org.mitk.gui.qt.cest/src/internal/QmitkCESTNormalizeView.cpp
@@ -1,110 +1,110 @@
/*============================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center (DKFZ)
All rights reserved.
Use of this source code is governed by a 3-clause BSD license that can be
found in the LICENSE file.
============================================================================*/
#include "QmitkCESTNormalizeView.h"
#include <QMessageBox>
#include "mitkWorkbenchUtil.h"
#include "mitkNodePredicateAnd.h"
#include "mitkNodePredicateDataProperty.h"
#include "mitkNodePredicateDataType.h"
#include "QmitkDataStorageComboBoxWithSelectNone.h"
#include <mitkImage.h>
#include "mitkCESTImageNormalizationFilter.h"
#include "mitkCustomTagParser.h"
#include "mitkCESTImageDetectionHelper.h"
const std::string QmitkCESTNormalizeView::VIEW_ID = "org.mitk.gui.qt.cest.normalize";
void QmitkCESTNormalizeView::SetFocus()
{
m_Controls.btnNormalize->setFocus();
}
void QmitkCESTNormalizeView::CreateQtPartControl(QWidget* parent)
{
m_Controls.setupUi(parent);
m_Controls.btnNormalize->setEnabled(false);
m_Controls.comboCESTImage->SetPredicate(this->m_IsCESTImagePredicate);
m_Controls.comboCESTImage->SetDataStorage(this->GetDataStorage());
connect(m_Controls.btnNormalize, SIGNAL(clicked()), this, SLOT(OnNormalizeButtonClicked()));
connect(m_Controls.comboCESTImage, SIGNAL(OnSelectionChanged(const mitk::DataNode *)), this, SLOT(UpdateGUIControls()));
UpdateGUIControls();
}
void QmitkCESTNormalizeView::UpdateGUIControls()
{
m_Controls.btnNormalize->setEnabled(m_Controls.comboCESTImage->GetSelectedNode().IsNotNull());
}
void QmitkCESTNormalizeView::OnNormalizeButtonClicked()
{
auto selectedImageNode = m_Controls.comboCESTImage->GetSelectedNode();
if (!selectedImageNode)
{
MITK_ERROR << "Invalid system state. CEST selection is invalid. Selected node is null_ptr.";
return;
}
auto selectedImage = dynamic_cast<mitk::Image*>(selectedImageNode->GetData());
if (!selectedImageNode)
{
MITK_ERROR << "Invalid system state. CEST selection is invalid. Selected node is not an image.";
return;
}
std::string offsetsStr = "";
- bool hasOffsets = selectedImage->GetPropertyList()->GetStringProperty(mitk::CustomTagParser::m_OffsetsPropertyName.c_str(), offsetsStr);
+ bool hasOffsets = selectedImage->GetPropertyList()->GetStringProperty(mitk::CEST_PROPERTY_NAME_OFFSETS().c_str(), offsetsStr);
if (!hasOffsets)
{
QMessageBox::information(nullptr, "CEST normalization", "Selected image was missing CEST offset information.");
return;
}
if (!mitk::IsNotNormalizedCESTImage(selectedImage))
{
QMessageBox::information(nullptr, "CEST normalization", "Selected image already seems to be normalized.");
return;
}
if (selectedImage->GetDimension() == 4)
{
auto normalizationFilter = mitk::CESTImageNormalizationFilter::New();
normalizationFilter->SetInput(selectedImage);
normalizationFilter->Update();
auto resultImage = normalizationFilter->GetOutput();
mitk::DataNode::Pointer dataNode = mitk::DataNode::New();
dataNode->SetData(resultImage);
std::string normalizedName = selectedImageNode->GetName() + "_normalized";
dataNode->SetName(normalizedName);
this->GetDataStorage()->Add(dataNode);
}
}
QmitkCESTNormalizeView::QmitkCESTNormalizeView()
{
auto isImage = mitk::NodePredicateDataType::New("Image");
this->m_IsCESTImagePredicate = mitk::NodePredicateAnd::New(isImage, mitk::CreateAnyCESTImageNodePredicate()).GetPointer();
}
diff --git a/Plugins/org.mitk.gui.qt.cest/src/internal/QmitkCESTStatisticsView.cpp b/Plugins/org.mitk.gui.qt.cest/src/internal/QmitkCESTStatisticsView.cpp
index 8c7f134f54..df0ace9532 100644
--- a/Plugins/org.mitk.gui.qt.cest/src/internal/QmitkCESTStatisticsView.cpp
+++ b/Plugins/org.mitk.gui.qt.cest/src/internal/QmitkCESTStatisticsView.cpp
@@ -1,816 +1,816 @@
/*============================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center (DKFZ)
All rights reserved.
Use of this source code is governed by a 3-clause BSD license that can be
found in the LICENSE file.
============================================================================*/
// itk
#include "itksys/SystemTools.hxx"
#include <itkImageRegionConstIterator.h>
#include <itkImageRegionIterator.h>
// Blueberry
#include <berryISelectionService.h>
#include <berryIWorkbenchWindow.h>
// Qmitk
#include "QmitkCESTStatisticsView.h"
// Qt
#include <QMessageBox>
#include <qclipboard.h>
// qwt
#include <qwt_scale_engine.h>
// mitk
#include <mitkCESTImageNormalizationFilter.h>
#include <mitkCustomTagParser.h>
#include <mitkITKImageImport.h>
#include <mitkImage.h>
#include <mitkImageAccessByItk.h>
#include <mitkImageCast.h>
#include <mitkImageStatisticsContainerNodeHelper.h>
#include <mitkLocaleSwitch.h>
#include <mitkSliceNavigationController.h>
#include <mitkStatisticsToImageRelationRule.h>
#include <mitkStatisticsToMaskRelationRule.h>
#include <mitkTemporoSpatialStringProperty.h>
#include <mitkTimeGeometry.h>
// boost
#include <boost/algorithm/string.hpp>
#include <boost/tokenizer.hpp>
// stl
#include <algorithm>
#include <iostream>
#include <iterator>
#include <sstream>
#include <string>
#include <vector>
namespace
{
template <typename T, typename Compare>
void GetSortPermutation(std::vector<unsigned> &out,
const std::vector<T> &determiningVector,
Compare compare = std::less<T>())
{
out.resize(determiningVector.size());
std::iota(out.begin(), out.end(), 0);
std::sort(out.begin(), out.end(), [&](unsigned i, unsigned j) {
return compare(determiningVector[i], determiningVector[j]);
});
}
template <typename T>
void ApplyPermutation(const std::vector<unsigned> &order, std::vector<T> &vectorToSort)
{
assert(order.size() == vectorToSort.size());
std::vector<T> tempVector(vectorToSort.size());
for (unsigned i = 0; i < vectorToSort.size(); i++)
{
tempVector[i] = vectorToSort[order[i]];
}
vectorToSort = tempVector;
}
template <typename T, typename... S>
void ApplyPermutation(const std::vector<unsigned> &order,
std::vector<T> ¤tVector,
std::vector<S> &... remainingVectors)
{
ApplyPermutation(order, currentVector);
ApplyPermutation(order, remainingVectors...);
}
template <typename T, typename Compare, typename... SS>
void SortVectors(const std::vector<T> &orderDeterminingVector,
Compare comparison,
std::vector<SS> &... vectorsToBeSorted)
{
std::vector<unsigned> order;
GetSortPermutation(order, orderDeterminingVector, comparison);
ApplyPermutation(order, vectorsToBeSorted...);
}
} // namespace
const std::string QmitkCESTStatisticsView::VIEW_ID = "org.mitk.views.ceststatistics";
QmitkCESTStatisticsView::QmitkCESTStatisticsView(QObject * /*parent*/, const char * /*name*/)
{
this->m_CalculatorJob = new QmitkImageStatisticsCalculationJob();
m_currentSelectedPosition.Fill(0.0);
m_currentSelectedTimeStep = 0;
m_CrosshairPointSet = mitk::PointSet::New();
}
QmitkCESTStatisticsView::~QmitkCESTStatisticsView()
{
while (this->m_CalculatorJob->isRunning()) // wait until thread has finished
{
itksys::SystemTools::Delay(100);
}
delete this->m_CalculatorJob;
}
void QmitkCESTStatisticsView::SetFocus()
{
m_Controls.threeDimToFourDimPushButton->setFocus();
}
void QmitkCESTStatisticsView::CreateQtPartControl(QWidget *parent)
{
// create GUI widgets from the Qt Designer's .ui file
m_Controls.setupUi(parent);
connect(
m_Controls.threeDimToFourDimPushButton, SIGNAL(clicked()), this, SLOT(OnThreeDimToFourDimPushButtonClicked()));
connect((QObject *)this->m_CalculatorJob,
SIGNAL(finished()),
this,
SLOT(OnThreadedStatisticsCalculationEnds()),
Qt::QueuedConnection);
connect((QObject *)(this->m_Controls.fixedRangeCheckBox),
SIGNAL(toggled(bool)),
(QObject *)this,
SLOT(OnFixedRangeCheckBoxToggled(bool)));
connect((QObject *)(this->m_Controls.fixedRangeLowerDoubleSpinBox),
SIGNAL(editingFinished()),
(QObject *)this,
SLOT(OnFixedRangeDoubleSpinBoxChanged()));
connect((QObject *)(this->m_Controls.fixedRangeUpperDoubleSpinBox),
SIGNAL(editingFinished()),
(QObject *)this,
SLOT(OnFixedRangeDoubleSpinBoxChanged()));
m_Controls.threeDimToFourDimPushButton->setEnabled(false);
m_Controls.widget_statistics->SetDataStorage(this->GetDataStorage());
this->m_SliceChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart());
connect(&m_SliceChangeListener, SIGNAL(SliceChanged()), this, SLOT(OnSliceChanged()));
}
void QmitkCESTStatisticsView::RenderWindowPartActivated(mitk::IRenderWindowPart *renderWindowPart)
{
this->m_SliceChangeListener.RenderWindowPartActivated(renderWindowPart);
}
void QmitkCESTStatisticsView::RenderWindowPartDeactivated(mitk::IRenderWindowPart *renderWindowPart)
{
this->m_SliceChangeListener.RenderWindowPartDeactivated(renderWindowPart);
}
void QmitkCESTStatisticsView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*source*/,
const QList<mitk::DataNode::Pointer> &nodes)
{
if (nodes.empty())
{
std::stringstream message;
message << "<font color='red'>Please select an image.</font>";
m_Controls.labelWarning->setText(message.str().c_str());
m_Controls.labelWarning->show();
this->Clear();
return;
}
// iterate all selected objects
bool atLeastOneWasCESTImage = false;
foreach (mitk::DataNode::Pointer node, nodes)
{
if (node.IsNull())
{
continue;
}
if (dynamic_cast<mitk::Image *>(node->GetData()) != nullptr)
{
m_Controls.labelWarning->setVisible(false);
bool zSpectrumSet = SetZSpectrum(dynamic_cast<mitk::StringProperty *>(
- node->GetData()->GetProperty(mitk::CustomTagParser::m_OffsetsPropertyName.c_str()).GetPointer()));
+ node->GetData()->GetProperty(mitk::CEST_PROPERTY_NAME_OFFSETS().c_str()).GetPointer()));
atLeastOneWasCESTImage = atLeastOneWasCESTImage || zSpectrumSet;
if (zSpectrumSet)
{
m_ZImage = dynamic_cast<mitk::Image *>(node->GetData());
m_Controls.widget_statistics->SetImageNodes({node.GetPointer()});
}
else
{
m_MaskImage = dynamic_cast<mitk::Image *>(node->GetData());
m_Controls.widget_statistics->SetMaskNodes({node.GetPointer()});
}
}
if (dynamic_cast<mitk::PlanarFigure *>(node->GetData()) != nullptr)
{
m_MaskPlanarFigure = dynamic_cast<mitk::PlanarFigure *>(node->GetData());
m_Controls.widget_statistics->SetMaskNodes({node.GetPointer()});
}
if (dynamic_cast<mitk::PointSet *>(node->GetData()) != nullptr)
{
m_PointSet = dynamic_cast<mitk::PointSet *>(node->GetData());
}
}
// We only want to offer normalization or timestep copying if one object is selected
if (nodes.size() == 1)
{
if (dynamic_cast<mitk::Image *>(nodes.front()->GetData()))
{
m_Controls.threeDimToFourDimPushButton->setDisabled(atLeastOneWasCESTImage);
}
else
{
m_Controls.threeDimToFourDimPushButton->setEnabled(false);
std::stringstream message;
message << "<font color='red'>The selected node is not an image.</font>";
m_Controls.labelWarning->setText(message.str().c_str());
m_Controls.labelWarning->show();
}
this->Clear();
return;
}
// we always need a mask, either image or planar figure as well as an image for further processing
if (nodes.size() != 2)
{
this->Clear();
return;
}
m_Controls.threeDimToFourDimPushButton->setEnabled(false);
if (!atLeastOneWasCESTImage)
{
std::stringstream message;
message << "<font color='red'>None of the selected data nodes contains required CEST meta information</font>";
m_Controls.labelWarning->setText(message.str().c_str());
m_Controls.labelWarning->show();
this->Clear();
return;
}
bool bothAreImages = (m_ZImage.GetPointer() != nullptr) && (m_MaskImage.GetPointer() != nullptr);
if (bothAreImages)
{
bool geometriesMatch =
mitk::Equal(*(m_ZImage->GetTimeGeometry()), *(m_MaskImage->GetTimeGeometry()), mitk::eps, false);
if (!geometriesMatch)
{
std::stringstream message;
message << "<font color='red'>The selected images have different geometries.</font>";
m_Controls.labelWarning->setText(message.str().c_str());
m_Controls.labelWarning->show();
this->Clear();
return;
}
}
if (!this->DataSanityCheck())
{
this->Clear();
return;
}
if (m_PointSet.IsNull())
{
// initialize thread and trigger it
this->m_CalculatorJob->SetIgnoreZeroValueVoxel(false);
this->m_CalculatorJob->Initialize(m_ZImage.GetPointer(), m_MaskImage.GetPointer(), m_MaskPlanarFigure.GetPointer());
std::stringstream message;
message << "<font color='red'>Calculating statistics...</font>";
m_Controls.labelWarning->setText(message.str().c_str());
m_Controls.labelWarning->show();
try
{
// Compute statistics
this->m_CalculatorJob->start();
}
catch (const mitk::Exception &e)
{
std::stringstream message;
message << "<font color='red'>" << e.GetDescription() << "</font>";
m_Controls.labelWarning->setText(message.str().c_str());
m_Controls.labelWarning->show();
}
catch (const std::runtime_error &e)
{
// In case of exception, print error message on GUI
std::stringstream message;
message << "<font color='red'>" << e.what() << "</font>";
m_Controls.labelWarning->setText(message.str().c_str());
m_Controls.labelWarning->show();
}
catch (const std::exception &e)
{
MITK_ERROR << "Caught exception: " << e.what();
// In case of exception, print error message on GUI
std::stringstream message;
message << "<font color='red'>Error! Unequal Dimensions of Image and Segmentation. No recompute possible </font>";
m_Controls.labelWarning->setText(message.str().c_str());
m_Controls.labelWarning->show();
}
while (this->m_CalculatorJob->isRunning()) // wait until thread has finished
{
itksys::SystemTools::Delay(100);
}
}
if (m_PointSet.IsNotNull())
{
if (m_ZImage->GetDimension() == 4)
{
AccessFixedDimensionByItk(m_ZImage, PlotPointSet, 4);
}
else
{
MITK_WARN << "Expecting a 4D image.";
}
}
}
void QmitkCESTStatisticsView::OnThreadedStatisticsCalculationEnds()
{
this->m_Controls.m_DataViewWidget->SetAxisTitle(QwtPlot::Axis::xBottom, "delta w");
this->m_Controls.m_DataViewWidget->SetAxisTitle(QwtPlot::Axis::yLeft, "z");
if (this->m_CalculatorJob->GetStatisticsUpdateSuccessFlag())
{
auto statistics = this->m_CalculatorJob->GetStatisticsData();
std::string statisticsNodeName = "CEST_statistics";
auto statisticsNode = mitk::CreateImageStatisticsNode(statistics, statisticsNodeName);
auto imageRule = mitk::StatisticsToImageRelationRule::New();
imageRule->Connect(statistics, m_CalculatorJob->GetStatisticsImage());
if (m_CalculatorJob->GetMaskImage())
{
auto maskRule = mitk::StatisticsToMaskRelationRule::New();
maskRule->Connect(statistics, m_CalculatorJob->GetMaskImage());
}
else if (m_CalculatorJob->GetPlanarFigure())
{
auto planarFigureRule = mitk::StatisticsToMaskRelationRule::New();
planarFigureRule->Connect(statistics, m_CalculatorJob->GetPlanarFigure());
}
this->GetDataStorage()->Add(statisticsNode);
QmitkPlotWidget::DataVector::size_type numberOfSpectra = this->m_zSpectrum.size();
QmitkPlotWidget::DataVector means(numberOfSpectra);
QmitkPlotWidget::DataVector stdevs(numberOfSpectra);
for (unsigned int index = 0; index < numberOfSpectra; ++index)
{
means[index] =
statistics->GetStatisticsForTimeStep(index).GetValueConverted<mitk::ImageStatisticsContainer::RealType>(
mitk::ImageStatisticsConstants::MEAN());
stdevs[index] =
statistics->GetStatisticsForTimeStep(index).GetValueConverted<mitk::ImageStatisticsContainer::RealType>(
mitk::ImageStatisticsConstants::STANDARDDEVIATION());
}
QmitkPlotWidget::DataVector xValues = this->m_zSpectrum;
RemoveMZeros(xValues, means, stdevs);
::SortVectors(xValues, std::less<qreal>(), xValues, means, stdevs);
unsigned int curveId = this->m_Controls.m_DataViewWidget->InsertCurve("Spectrum");
this->m_Controls.m_DataViewWidget->SetCurveData(curveId, xValues, means, stdevs, stdevs);
this->m_Controls.m_DataViewWidget->SetErrorPen(curveId, QPen(Qt::blue));
QwtSymbol *blueSymbol = new QwtSymbol(QwtSymbol::Rect, QColor(Qt::blue), QColor(Qt::blue), QSize(8, 8));
this->m_Controls.m_DataViewWidget->SetCurveSymbol(curveId, blueSymbol);
this->m_Controls.m_DataViewWidget->SetLegendAttribute(curveId, QwtPlotCurve::LegendShowSymbol);
QwtLegend *legend = new QwtLegend();
legend->setFrameShape(QFrame::Box);
legend->setFrameShadow(QFrame::Sunken);
legend->setLineWidth(1);
this->m_Controls.m_DataViewWidget->SetLegend(legend, QwtPlot::BottomLegend);
m_Controls.m_DataViewWidget->GetPlot()
->axisScaleEngine(QwtPlot::Axis::xBottom)
->setAttributes(QwtScaleEngine::Inverted);
this->m_Controls.m_DataViewWidget->Replot();
m_Controls.labelWarning->setVisible(false);
m_Controls.m_StatisticsGroupBox->setEnabled(true);
m_Controls.m_StatisticsGroupBox->setEnabled(true);
if (this->m_Controls.fixedRangeCheckBox->isChecked())
{
this->m_Controls.m_DataViewWidget->GetPlot()->setAxisAutoScale(2, false);
this->m_Controls.m_DataViewWidget->GetPlot()->setAxisScale(
2,
this->m_Controls.fixedRangeLowerDoubleSpinBox->value(),
this->m_Controls.fixedRangeUpperDoubleSpinBox->value());
}
else
{
this->m_Controls.m_DataViewWidget->GetPlot()->setAxisAutoScale(2, true);
}
}
else
{
m_Controls.labelWarning->setText(m_CalculatorJob->GetLastErrorMessage().c_str());
m_Controls.labelWarning->setVisible(true);
this->Clear();
}
}
void QmitkCESTStatisticsView::OnFixedRangeDoubleSpinBoxChanged()
{
if (this->m_Controls.fixedRangeCheckBox->isChecked())
{
this->m_Controls.m_DataViewWidget->GetPlot()->setAxisAutoScale(2, false);
this->m_Controls.m_DataViewWidget->GetPlot()->setAxisScale(2,
this->m_Controls.fixedRangeLowerDoubleSpinBox->value(),
this->m_Controls.fixedRangeUpperDoubleSpinBox->value());
}
this->m_Controls.m_DataViewWidget->Replot();
}
template <typename TPixel, unsigned int VImageDimension>
void QmitkCESTStatisticsView::PlotPointSet(itk::Image<TPixel, VImageDimension> *image)
{
this->m_Controls.m_DataViewWidget->SetAxisTitle(QwtPlot::Axis::xBottom, "delta w");
this->m_Controls.m_DataViewWidget->SetAxisTitle(QwtPlot::Axis::yLeft, "z");
QmitkPlotWidget::DataVector::size_type numberOfSpectra = this->m_zSpectrum.size();
mitk::PointSet::Pointer internalPointset;
if (m_PointSet.IsNotNull())
{
internalPointset = m_PointSet;
}
else
{
internalPointset = m_CrosshairPointSet;
}
if (internalPointset.IsNull())
{
return;
}
if (!this->DataSanityCheck())
{
m_Controls.labelWarning->setText("Data can not be plotted, internally inconsistent.");
m_Controls.labelWarning->show();
return;
}
auto maxIndex = internalPointset->GetMaxId().Index();
for (std::size_t number = 0; number < maxIndex + 1; ++number)
{
mitk::PointSet::PointType point;
if (!internalPointset->GetPointIfExists(number, &point))
{
continue;
}
if (!this->m_ZImage->GetGeometry()->IsInside(point))
{
continue;
}
itk::Index<3> itkIndex;
this->m_ZImage->GetGeometry()->WorldToIndex(point, itkIndex);
itk::Index<VImageDimension> itkIndexTime;
itkIndexTime[0] = itkIndex[0];
itkIndexTime[1] = itkIndex[1];
itkIndexTime[2] = itkIndex[2];
QmitkPlotWidget::DataVector values(numberOfSpectra);
for (std::size_t step = 0; step < numberOfSpectra; ++step)
{
if (VImageDimension == 4)
{
itkIndexTime[3] = step;
}
values[step] = image->GetPixel(itkIndexTime);
}
std::stringstream name;
name << "Point " << number;
// Qcolor enums go from 0 to 19, but 19 is transparent and 0,1 are for bitmaps
// 3 is white and thus not visible
QColor color(static_cast<Qt::GlobalColor>(number % 17 + 4));
QmitkPlotWidget::DataVector xValues = this->m_zSpectrum;
RemoveMZeros(xValues, values);
::SortVectors(xValues, std::less<qreal>(), xValues, values);
unsigned int curveId = this->m_Controls.m_DataViewWidget->InsertCurve(name.str().c_str());
this->m_Controls.m_DataViewWidget->SetCurveData(curveId, xValues, values);
this->m_Controls.m_DataViewWidget->SetCurvePen(curveId, QPen(color));
QwtSymbol *symbol = new QwtSymbol(QwtSymbol::Rect, color, color, QSize(8, 8));
this->m_Controls.m_DataViewWidget->SetCurveSymbol(curveId, symbol);
this->m_Controls.m_DataViewWidget->SetLegendAttribute(curveId, QwtPlotCurve::LegendShowSymbol);
}
if (this->m_Controls.fixedRangeCheckBox->isChecked())
{
this->m_Controls.m_DataViewWidget->GetPlot()->setAxisAutoScale(2, false);
this->m_Controls.m_DataViewWidget->GetPlot()->setAxisScale(2,
this->m_Controls.fixedRangeLowerDoubleSpinBox->value(),
this->m_Controls.fixedRangeUpperDoubleSpinBox->value());
}
else
{
this->m_Controls.m_DataViewWidget->GetPlot()->setAxisAutoScale(2, true);
}
QwtLegend *legend = new QwtLegend();
legend->setFrameShape(QFrame::Box);
legend->setFrameShadow(QFrame::Sunken);
legend->setLineWidth(1);
this->m_Controls.m_DataViewWidget->SetLegend(legend, QwtPlot::BottomLegend);
m_Controls.m_DataViewWidget->GetPlot()
->axisScaleEngine(QwtPlot::Axis::xBottom)
->setAttributes(QwtScaleEngine::Inverted);
this->m_Controls.m_DataViewWidget->Replot();
m_Controls.labelWarning->setVisible(false);
}
void QmitkCESTStatisticsView::OnFixedRangeCheckBoxToggled(bool state)
{
this->m_Controls.fixedRangeLowerDoubleSpinBox->setEnabled(state);
this->m_Controls.fixedRangeUpperDoubleSpinBox->setEnabled(state);
}
void QmitkCESTStatisticsView::RemoveMZeros(QmitkPlotWidget::DataVector &xValues, QmitkPlotWidget::DataVector &yValues)
{
QmitkPlotWidget::DataVector tempX;
QmitkPlotWidget::DataVector tempY;
for (std::size_t index = 0; index < xValues.size(); ++index)
{
if ((xValues.at(index) < -299) || (xValues.at(index)) > 299)
{
// do not include
}
else
{
tempX.push_back(xValues.at(index));
tempY.push_back(yValues.at(index));
}
}
xValues = tempX;
yValues = tempY;
}
void QmitkCESTStatisticsView::RemoveMZeros(QmitkPlotWidget::DataVector &xValues,
QmitkPlotWidget::DataVector &yValues,
QmitkPlotWidget::DataVector &stdDevs)
{
QmitkPlotWidget::DataVector tempX;
QmitkPlotWidget::DataVector tempY;
QmitkPlotWidget::DataVector tempDevs;
for (std::size_t index = 0; index < xValues.size(); ++index)
{
if ((xValues.at(index) < -299) || (xValues.at(index)) > 299)
{
// do not include
}
else
{
tempX.push_back(xValues.at(index));
tempY.push_back(yValues.at(index));
tempDevs.push_back(stdDevs.at(index));
}
}
xValues = tempX;
yValues = tempY;
stdDevs = tempDevs;
}
void QmitkCESTStatisticsView::OnThreeDimToFourDimPushButtonClicked()
{
QList<mitk::DataNode::Pointer> nodes = this->GetDataManagerSelection();
if (nodes.empty())
return;
mitk::DataNode *node = nodes.front();
if (!node)
{
// Nothing selected. Inform the user and return
QMessageBox::information(nullptr, "CEST View", "Please load and select an image before starting image processing.");
return;
}
// here we have a valid mitk::DataNode
// a node itself is not very useful, we need its data item (the image)
mitk::BaseData *data = node->GetData();
if (data)
{
// test if this data item is an image or not (could also be a surface or something totally different)
mitk::Image *image = dynamic_cast<mitk::Image *>(data);
if (image)
{
if (image->GetDimension() == 4)
{
AccessFixedDimensionByItk(image, CopyTimesteps, 4);
}
this->Clear();
}
}
}
template <typename TPixel, unsigned int VImageDimension>
void QmitkCESTStatisticsView::CopyTimesteps(itk::Image<TPixel, VImageDimension> *image)
{
typedef itk::Image<TPixel, VImageDimension> ImageType;
// typedef itk::PasteImageFilter<ImageType, ImageType> PasteImageFilterType;
unsigned int numberOfTimesteps = image->GetLargestPossibleRegion().GetSize(3);
typename ImageType::RegionType sourceRegion = image->GetLargestPossibleRegion();
sourceRegion.SetSize(3, 1);
typename ImageType::RegionType targetRegion = image->GetLargestPossibleRegion();
targetRegion.SetSize(3, 1);
for (unsigned int timestep = 1; timestep < numberOfTimesteps; ++timestep)
{
targetRegion.SetIndex(3, timestep);
itk::ImageRegionConstIterator<ImageType> sourceIterator(image, sourceRegion);
itk::ImageRegionIterator<ImageType> targetIterator(image, targetRegion);
while (!sourceIterator.IsAtEnd())
{
targetIterator.Set(sourceIterator.Get());
++sourceIterator;
++targetIterator;
}
}
}
bool QmitkCESTStatisticsView::SetZSpectrum(mitk::StringProperty *zSpectrumProperty)
{
if (zSpectrumProperty == nullptr)
{
return false;
}
mitk::LocaleSwitch localeSwitch("C");
std::string zSpectrumString = zSpectrumProperty->GetValueAsString();
std::istringstream iss(zSpectrumString);
std::vector<std::string> zSpectra;
std::copy(
std::istream_iterator<std::string>(iss), std::istream_iterator<std::string>(), std::back_inserter(zSpectra));
m_zSpectrum.clear();
m_zSpectrum.resize(0);
for (auto const &spectrumString : zSpectra)
{
m_zSpectrum.push_back(std::stod(spectrumString));
}
return (m_zSpectrum.size() > 0);
}
bool QmitkCESTStatisticsView::DataSanityCheck()
{
QmitkPlotWidget::DataVector::size_type numberOfSpectra = m_zSpectrum.size();
// if we do not have a spectrum, the data can not be processed
if (numberOfSpectra == 0)
{
return false;
}
// if we do not have CEST image data, the data can not be processed
if (m_ZImage.IsNull())
{
return false;
}
// if the CEST image data and the meta information do not match, the data can not be processed
if (numberOfSpectra != m_ZImage->GetTimeSteps())
{
MITK_INFO << "CEST meta information and number of volumes does not match.";
return false;
}
// if we have neither a mask image, a point set nor a mask planar figure, we can not do statistics
// statistics on the whole image would not make sense
if (m_MaskImage.IsNull() && m_MaskPlanarFigure.IsNull() && m_PointSet.IsNull() && m_CrosshairPointSet->IsEmpty())
{
return false;
}
// if we have a mask image and a mask planar figure, we can not do statistics
// we do not know which one to use
if (m_MaskImage.IsNotNull() && m_MaskPlanarFigure.IsNotNull())
{
return false;
}
return true;
}
void QmitkCESTStatisticsView::Clear()
{
this->m_zSpectrum.clear();
this->m_zSpectrum.resize(0);
this->m_ZImage = nullptr;
this->m_MaskImage = nullptr;
this->m_MaskPlanarFigure = nullptr;
this->m_PointSet = nullptr;
this->m_Controls.m_DataViewWidget->Clear();
this->m_Controls.m_StatisticsGroupBox->setEnabled(false);
this->m_Controls.widget_statistics->SetImageNodes({});
this->m_Controls.widget_statistics->SetMaskNodes({});
}
void QmitkCESTStatisticsView::OnSliceChanged()
{
mitk::Point3D currentSelectedPosition = this->GetRenderWindowPart()->GetSelectedPosition(nullptr);
unsigned int currentSelectedTimeStep =
this->GetRenderWindowPart()->GetTimeNavigationController()->GetTime()->GetPos();
if (m_currentSelectedPosition != currentSelectedPosition || m_currentSelectedTimeStep != currentSelectedTimeStep)
//|| m_selectedNodeTime > m_currentPositionTime)
{
// the current position has been changed or the selected node has been changed since the last position validation ->
// check position
m_currentSelectedPosition = currentSelectedPosition;
m_currentSelectedTimeStep = currentSelectedTimeStep;
m_currentPositionTime.Modified();
m_CrosshairPointSet->Clear();
m_CrosshairPointSet->SetPoint(0, m_currentSelectedPosition);
QList<mitk::DataNode::Pointer> nodes = this->GetDataManagerSelection();
if (nodes.empty() || nodes.size() > 1)
return;
mitk::DataNode *node = nodes.front();
if (!node)
{
return;
}
if (dynamic_cast<mitk::Image *>(node->GetData()) != nullptr)
{
m_Controls.labelWarning->setVisible(false);
bool zSpectrumSet = SetZSpectrum(dynamic_cast<mitk::StringProperty *>(
- node->GetData()->GetProperty(mitk::CustomTagParser::m_OffsetsPropertyName.c_str()).GetPointer()));
+ node->GetData()->GetProperty(mitk::CEST_PROPERTY_NAME_OFFSETS().c_str()).GetPointer()));
if (zSpectrumSet)
{
m_ZImage = dynamic_cast<mitk::Image *>(node->GetData());
}
else
{
return;
}
}
else
{
return;
}
this->m_Controls.m_DataViewWidget->Clear();
AccessFixedDimensionByItk(m_ZImage, PlotPointSet, 4);
}
}