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 + +#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 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 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 +#include #include #include #include #include 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 ExtractOffsets(const mitk::Image* image) -{ - std::vector result; - - if (image) - { - std::string offsets = ""; - std::vector 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 void mitk::CESTImageNormalizationFilter::NormalizeTimeSteps(const itk::Image* image) { typedef itk::Image ImageType; typedef itk::Image OutputImageType; - auto offsets = ExtractOffsets(this->GetInput()); + auto offsets = ExtractCESTOffset(this->GetInput()); // determine normalization images std::vector 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 lowerMZeroIterator(image, lowerMZeroRegion); itk::ImageRegionConstIterator upperMZeroIterator(image, upperMZeroRegion); itk::ImageRegionConstIterator sourceIterator(image, sourceRegion); itk::ImageRegionIterator 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(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 +#include + + +std::vector mitk::ExtractCESTT1Time(const BaseData* image) +{ + std::vector 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::ExtractCESTOffset(const BaseData* image) +{ + std::vector 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 #include "mitkWorkbenchUtil.h" #include "mitkNodePredicateAnd.h" #include "mitkNodePredicateDataProperty.h" #include "mitkNodePredicateDataType.h" #include "QmitkDataStorageComboBoxWithSelectNone.h" #include #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(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 #include // Blueberry #include #include // Qmitk #include "QmitkCESTStatisticsView.h" // Qt #include #include // qwt #include // mitk #include #include #include #include #include #include #include #include #include #include #include #include #include // boost #include #include // stl #include #include #include #include #include #include namespace { template void GetSortPermutation(std::vector &out, const std::vector &determiningVector, Compare compare = std::less()) { 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 void ApplyPermutation(const std::vector &order, std::vector &vectorToSort) { assert(order.size() == vectorToSort.size()); std::vector tempVector(vectorToSort.size()); for (unsigned i = 0; i < vectorToSort.size(); i++) { tempVector[i] = vectorToSort[order[i]]; } vectorToSort = tempVector; } template void ApplyPermutation(const std::vector &order, std::vector ¤tVector, std::vector &... remainingVectors) { ApplyPermutation(order, currentVector); ApplyPermutation(order, remainingVectors...); } template void SortVectors(const std::vector &orderDeterminingVector, Compare comparison, std::vector &... vectorsToBeSorted) { std::vector 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 &nodes) { if (nodes.empty()) { std::stringstream message; message << "Please select an image."; 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(node->GetData()) != nullptr) { m_Controls.labelWarning->setVisible(false); bool zSpectrumSet = SetZSpectrum(dynamic_cast( - 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(node->GetData()); m_Controls.widget_statistics->SetImageNodes({node.GetPointer()}); } else { m_MaskImage = dynamic_cast(node->GetData()); m_Controls.widget_statistics->SetMaskNodes({node.GetPointer()}); } } if (dynamic_cast(node->GetData()) != nullptr) { m_MaskPlanarFigure = dynamic_cast(node->GetData()); m_Controls.widget_statistics->SetMaskNodes({node.GetPointer()}); } if (dynamic_cast(node->GetData()) != nullptr) { m_PointSet = dynamic_cast(node->GetData()); } } // We only want to offer normalization or timestep copying if one object is selected if (nodes.size() == 1) { if (dynamic_cast(nodes.front()->GetData())) { m_Controls.threeDimToFourDimPushButton->setDisabled(atLeastOneWasCESTImage); } else { m_Controls.threeDimToFourDimPushButton->setEnabled(false); std::stringstream message; message << "The selected node is not an image."; 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 << "None of the selected data nodes contains required CEST meta information"; 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 << "The selected images have different geometries."; 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 << "Calculating statistics..."; 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 << "" << e.GetDescription() << ""; 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 << "" << e.what() << ""; 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 << "Error! Unequal Dimensions of Image and Segmentation. No recompute possible "; 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::ImageStatisticsConstants::MEAN()); stdevs[index] = statistics->GetStatisticsForTimeStep(index).GetValueConverted( mitk::ImageStatisticsConstants::STANDARDDEVIATION()); } QmitkPlotWidget::DataVector xValues = this->m_zSpectrum; RemoveMZeros(xValues, means, stdevs); ::SortVectors(xValues, std::less(), 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 void QmitkCESTStatisticsView::PlotPointSet(itk::Image *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 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(number % 17 + 4)); QmitkPlotWidget::DataVector xValues = this->m_zSpectrum; RemoveMZeros(xValues, values); ::SortVectors(xValues, std::less(), 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 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(data); if (image) { if (image->GetDimension() == 4) { AccessFixedDimensionByItk(image, CopyTimesteps, 4); } this->Clear(); } } } template void QmitkCESTStatisticsView::CopyTimesteps(itk::Image *image) { typedef itk::Image ImageType; // typedef itk::PasteImageFilter 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 sourceIterator(image, sourceRegion); itk::ImageRegionIterator 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 zSpectra; std::copy( std::istream_iterator(iss), std::istream_iterator(), 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 nodes = this->GetDataManagerSelection(); if (nodes.empty() || nodes.size() > 1) return; mitk::DataNode *node = nodes.front(); if (!node) { return; } if (dynamic_cast(node->GetData()) != nullptr) { m_Controls.labelWarning->setVisible(false); bool zSpectrumSet = SetZSpectrum(dynamic_cast( - 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(node->GetData()); } else { return; } } else { return; } this->m_Controls.m_DataViewWidget->Clear(); AccessFixedDimensionByItk(m_ZImage, PlotPointSet, 4); } }