diff --git a/Modules/ImageStatistics/files.cmake b/Modules/ImageStatistics/files.cmake index 7b16237156..11ed3564bc 100644 --- a/Modules/ImageStatistics/files.cmake +++ b/Modules/ImageStatistics/files.cmake @@ -1,6 +1,7 @@ set(CPP_FILES mitkImageStatisticsCalculator.cpp mitkPointSetStatisticsCalculator.cpp mitkPointSetDifferenceStatisticsCalculator.cpp + mitkIntensityProfile.cpp ) diff --git a/Modules/ImageStatistics/mitkIntensityProfile.cpp b/Modules/ImageStatistics/mitkIntensityProfile.cpp new file mode 100644 index 0000000000..9199a3bc78 --- /dev/null +++ b/Modules/ImageStatistics/mitkIntensityProfile.cpp @@ -0,0 +1,311 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include +#include +#include +#include +#include +#include +#include +#include "mitkIntensityProfile.h" + +using namespace mitk; + +template +static void ReadPixel(const PixelType&, Image::Pointer image, const Index3D& index, ScalarType* returnValue) +{ + switch (image->GetDimension()) + { + case 2: + { + ImagePixelReadAccessor readAccess(image, image->GetSliceData(0)); + *returnValue = readAccess.GetPixelByIndex(reinterpret_cast&>(index)); + break; + } + + case 3: + { + ImagePixelReadAccessor readAccess(image, image->GetVolumeData(0)); + *returnValue = readAccess.GetPixelByIndex(index); + break; + } + + default: + *returnValue = 0; + break; + } +} + +static IntensityProfile::Pointer ComputeIntensityProfile(Image::Pointer image, itk::PolyLineParametricPath<3>::Pointer path) +{ + IntensityProfile::Pointer intensityProfile = IntensityProfile::New(); + itk::PolyLineParametricPath<3>::InputType input = path->StartOfInput(); + Geometry3D* imageGeometry = image->GetGeometry(); + const PixelType pixelType = image->GetPixelType(); + + IntensityProfile::MeasurementVectorType measurementVector; + itk::PolyLineParametricPath<3>::OffsetType offset; + Point3D worldPoint; + Index3D index; + + do + { + imageGeometry->IndexToWorld(path->Evaluate(input), worldPoint); + imageGeometry->WorldToIndex(worldPoint, index); + + mitkPixelTypeMultiplex3(ReadPixel, pixelType, image, index, measurementVector.GetDataPointer()); + intensityProfile->PushBack(measurementVector); + + offset = path->IncrementInput(input); + } while ((offset[0] | offset[1] | offset[2]) != 0); + + return intensityProfile; +} + +template +static typename itk::InterpolateImageFunction::Pointer CreateInterpolateImageFunction(InterpolateImageFunction::Enum interpolator) +{ + switch (interpolator) + { + case InterpolateImageFunction::NearestNeighbor: + return itk::NearestNeighborInterpolateImageFunction::New().GetPointer(); + + case InterpolateImageFunction::Linear: + return itk::LinearInterpolateImageFunction::New().GetPointer(); + + case InterpolateImageFunction::WindowedSinc_Blackman_3: + return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); + + case InterpolateImageFunction::WindowedSinc_Blackman_4: + return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); + + case InterpolateImageFunction::WindowedSinc_Blackman_5: + return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); + + case InterpolateImageFunction::WindowedSinc_Cosine_3: + return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); + + case InterpolateImageFunction::WindowedSinc_Cosine_4: + return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); + + case InterpolateImageFunction::WindowedSinc_Cosine_5: + return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); + + case InterpolateImageFunction::WindowedSinc_Hamming_3: + return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); + + case InterpolateImageFunction::WindowedSinc_Hamming_4: + return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); + + case InterpolateImageFunction::WindowedSinc_Hamming_5: + return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); + + case InterpolateImageFunction::WindowedSinc_Lanczos_3: + return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); + + case InterpolateImageFunction::WindowedSinc_Lanczos_4: + return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); + + case InterpolateImageFunction::WindowedSinc_Lanczos_5: + return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); + + case InterpolateImageFunction::WindowedSinc_Welch_3: + return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); + + case InterpolateImageFunction::WindowedSinc_Welch_4: + return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); + + case InterpolateImageFunction::WindowedSinc_Welch_5: + return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); + + default: + return itk::NearestNeighborInterpolateImageFunction::New().GetPointer(); + } +} + +template +static void ComputeIntensityProfile(itk::Image* image, itk::PolyLineParametricPath<3>::Pointer path, unsigned int numSamples, InterpolateImageFunction::Enum interpolator, IntensityProfile::Pointer intensityProfile) +{ + typename itk::InterpolateImageFunction >::Pointer interpolateImageFunction = CreateInterpolateImageFunction >(interpolator); + interpolateImageFunction->SetInputImage(image); + + const itk::PolyLineParametricPath<3>::InputType startOfInput = path->StartOfInput(); + const itk::PolyLineParametricPath<3>::InputType delta = 1.0 / (numSamples - 1); + + IntensityProfile::MeasurementVectorType measurementVector; + + for (unsigned int i = 0; i < numSamples; ++i) + { + measurementVector[0] = interpolateImageFunction->EvaluateAtContinuousIndex(path->Evaluate(startOfInput + i * delta)); + intensityProfile->PushBack(measurementVector); + } +} + +static IntensityProfile::Pointer ComputeIntensityProfile(Image::Pointer image, itk::PolyLineParametricPath<3>::Pointer path, unsigned int numSamples, InterpolateImageFunction::Enum interpolator) +{ + IntensityProfile::Pointer intensityProfile = IntensityProfile::New(); + AccessFixedDimensionByItk_n(image, ComputeIntensityProfile, 3, (path, numSamples, interpolator, intensityProfile)); + return intensityProfile; +} + +class AddPointToPath +{ +public: + AddPointToPath(const Geometry2D* planarFigureGeometry, const Geometry3D* imageGeometry, itk::PolyLineParametricPath<3>::Pointer path) + : m_PlanarFigureGeometry(planarFigureGeometry), + m_ImageGeometry(imageGeometry), + m_Path(path) + { + } + + void operator()(const PlanarFigure::PolyLineElement& polyLineElement) + { + m_PlanarFigureGeometry->Map(polyLineElement.Point, m_WorldPoint); + m_ImageGeometry->WorldToIndex(m_WorldPoint, m_ContinuousIndexPoint); + m_Vertex.CastFrom(m_ContinuousIndexPoint); + m_Path->AddVertex(m_Vertex); + } + +private: + const Geometry2D* m_PlanarFigureGeometry; + const Geometry3D* m_ImageGeometry; + itk::PolyLineParametricPath<3>::Pointer m_Path; + + Point3D m_WorldPoint; + Point3D m_ContinuousIndexPoint; + itk::PolyLineParametricPath<3>::ContinuousIndexType m_Vertex; +}; + +static itk::PolyLineParametricPath<3>::Pointer CreatePathFromPlanarFigure(Geometry3D* imageGeometry, PlanarFigure* planarFigure) +{ + itk::PolyLineParametricPath<3>::Pointer path = itk::PolyLineParametricPath<3>::New(); + const PlanarFigure::PolyLineType polyLine = planarFigure->GetPolyLine(0); + + std::for_each(polyLine.begin(), polyLine.end(), + AddPointToPath(planarFigure->GetGeometry2D(), imageGeometry, path)); + + return path; +} + +IntensityProfile::Pointer mitk::ComputeIntensityProfile(Image::Pointer image, PlanarFigure::Pointer planarFigure) +{ + return ::ComputeIntensityProfile(image, CreatePathFromPlanarFigure(image->GetGeometry(), planarFigure)); +} + +IntensityProfile::Pointer mitk::ComputeIntensityProfile(Image::Pointer image, PlanarLine::Pointer planarLine, unsigned int numSamples, InterpolateImageFunction::Enum interpolator) +{ + return ::ComputeIntensityProfile(image, CreatePathFromPlanarFigure(image->GetGeometry(), planarLine.GetPointer()), numSamples, interpolator); +} + +IntensityProfile::InstanceIdentifier mitk::ComputeGlobalMaximum(IntensityProfile::Pointer intensityProfile) +{ + IntensityProfile::MeasurementType max = -vcl_numeric_limits::max(); + IntensityProfile::InstanceIdentifier maxIndex = 0; + + IntensityProfile::ConstIterator end = intensityProfile->End(); + IntensityProfile::MeasurementType measurement; + + for (IntensityProfile::ConstIterator it = intensityProfile->Begin(); it != end; ++it) + { + measurement = it.GetMeasurementVector()[0]; + + if (measurement > max) + { + max = measurement; + maxIndex = it.GetInstanceIdentifier(); + } + } + + return maxIndex; +} + +IntensityProfile::InstanceIdentifier mitk::ComputeGlobalMinimum(IntensityProfile::Pointer intensityProfile) +{ + IntensityProfile::MeasurementType min = vcl_numeric_limits::max(); + IntensityProfile::InstanceIdentifier minIndex = 0; + + IntensityProfile::ConstIterator end = intensityProfile->End(); + IntensityProfile::MeasurementType measurement; + + for (IntensityProfile::ConstIterator it = intensityProfile->Begin(); it != end; ++it) + { + measurement = it.GetMeasurementVector()[0]; + + if (measurement < min) + { + min = measurement; + minIndex = it.GetInstanceIdentifier(); + } + } + + return minIndex; +} + +IntensityProfile::InstanceIdentifier mitk::ComputeCenterOfMaximumArea(IntensityProfile::Pointer intensityProfile, IntensityProfile::InstanceIdentifier radius) +{ + const IntensityProfile::MeasurementType min = intensityProfile->GetMeasurementVector(ComputeGlobalMinimum(intensityProfile))[0]; + const IntensityProfile::InstanceIdentifier areaWidth = 1 + 2 * radius; + + IntensityProfile::MeasurementType maxArea = 0; + + for (IntensityProfile::InstanceIdentifier i = 0; i < areaWidth; ++i) + maxArea += intensityProfile->GetMeasurementVector(i)[0] - min; + + const IntensityProfile::InstanceIdentifier lastIndex = intensityProfile->Size() - areaWidth; + IntensityProfile::InstanceIdentifier centerOfMaxArea = radius; + IntensityProfile::MeasurementType area = maxArea; + + for (IntensityProfile::InstanceIdentifier i = 1; i <= lastIndex; ++i) + { + area += intensityProfile->GetMeasurementVector(i + areaWidth - 1)[0] - min; + area -= intensityProfile->GetMeasurementVector(i - 1)[0] - min; + + if (area > maxArea) + { + maxArea = area; + centerOfMaxArea = i + radius; // TODO: If multiple areas in the neighborhood have the same intensity chose the middle one instead of the first one. + } + } + + return centerOfMaxArea; +} + +std::vector mitk::CreateVectorFromIntensityProfile(IntensityProfile::Pointer intensityProfile) +{ + std::vector result; + result.reserve(intensityProfile->Size()); + + IntensityProfile::ConstIterator end = intensityProfile->End(); + + for (IntensityProfile::ConstIterator it = intensityProfile->Begin(); it != end; ++it) + result.push_back(it.GetMeasurementVector()[0]); + + return result; +} + +IntensityProfile::Pointer mitk::CreateIntensityProfileFromVector(const std::vector& vector) +{ + const IntensityProfile::InstanceIdentifier size = vector.size(); + + IntensityProfile::Pointer result = IntensityProfile::New(); + result->Resize(size); + + for (IntensityProfile::InstanceIdentifier i = 0; i < size; ++i) + result->SetMeasurement(i, 0, vector[i]); + + return result; +} diff --git a/Modules/ImageStatistics/mitkIntensityProfile.h b/Modules/ImageStatistics/mitkIntensityProfile.h new file mode 100644 index 0000000000..f4086acc9c --- /dev/null +++ b/Modules/ImageStatistics/mitkIntensityProfile.h @@ -0,0 +1,115 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#ifndef mitkIntensityProfile_h +#define mitkIntensityProfile_h + +#include +#include +#include +#include + +namespace mitk +{ + typedef itk::Statistics::ListSample::MeasurementVectorType> IntensityProfile; + + /** \brief Compute intensity profile of an image for each pixel along the first PolyLine of a given planar figure. + * + * \param[in] image A two or three-dimensional image which consists of single component pixels. + * \param[in] planarFigure A planar figure from which the first PolyLine is used to evaluate the intensity profile. + * + * \return The computed intensity profile. + */ + MitkImageStatistics_EXPORT IntensityProfile::Pointer ComputeIntensityProfile(Image::Pointer image, PlanarFigure::Pointer planarFigure); + + namespace InterpolateImageFunction + { + enum Enum + { + NearestNeighbor, + Linear, + WindowedSinc_Blackman_3, + WindowedSinc_Blackman_4, + WindowedSinc_Blackman_5, + WindowedSinc_Cosine_3, + WindowedSinc_Cosine_4, + WindowedSinc_Cosine_5, + WindowedSinc_Hamming_3, + WindowedSinc_Hamming_4, + WindowedSinc_Hamming_5, + WindowedSinc_Lanczos_3, + WindowedSinc_Lanczos_4, + WindowedSinc_Lanczos_5, + WindowedSinc_Welch_3, + WindowedSinc_Welch_4, + WindowedSinc_Welch_5 + }; + } + + /** \brief Compute intensity profile of an image for each sample along a planar line. + * + * \param[in] image A three-dimensional image which consists of single component pixels. + * \param[in] planarLine A planar line along which the intensity profile will be evaluated. + * \param[in] numSamples Number of samples along the planar line (must be at least 2). + * \param[in] interpolator Image interpolation function which is used to read each sample. + * + * \return The computed intensity profile. + */ + MitkImageStatistics_EXPORT IntensityProfile::Pointer ComputeIntensityProfile(Image::Pointer image, PlanarLine::Pointer planarLine, unsigned int numSamples, InterpolateImageFunction::Enum interpolator = InterpolateImageFunction::NearestNeighbor); + + /** \brief Compute global maximum of an intensity profile. + * + * \param[in] intensityProfile An intensity profile. + * + * \return Index of the global maximum. + */ + MitkImageStatistics_EXPORT IntensityProfile::InstanceIdentifier ComputeGlobalMaximum(IntensityProfile::Pointer intensityProfile); + + /** \brief Compute global minimum of an intensity profile. + * + * \param[in] intensityProfile An intensity profile. + * + * \return Index of the global minimum. + */ + MitkImageStatistics_EXPORT IntensityProfile::InstanceIdentifier ComputeGlobalMinimum(IntensityProfile::Pointer intensityProfile); + + /** \brief Compute center of maximum area under the curve of an intensity profile. + * + * \param[in] intensityProfile An intensity profile. + * \param[in] radius Radius of the area (width of area equals 1 + 2 * radius). + * + * \return Index of the maximum area center. + */ + MitkImageStatistics_EXPORT IntensityProfile::InstanceIdentifier ComputeCenterOfMaximumArea(IntensityProfile::Pointer intensityProfile, IntensityProfile::InstanceIdentifier radius); + + /** \brief Convert an intensity profile to a standard library vector. + * + * \param[in] intensityProfile An intensity profile. + * + * \return Standard library vector which contains the input intensity profile measurements. + */ + MitkImageStatistics_EXPORT std::vector CreateVectorFromIntensityProfile(IntensityProfile::Pointer intensityProfile); + + /** \brief Convert a standard library vector to an intensity profile. + * + * \param[in] vector An standard library vector which contains intensity profile measurements. + * + * \return An intensity profile. + */ + MitkImageStatistics_EXPORT IntensityProfile::Pointer CreateIntensityProfileFromVector(const std::vector& vector); +} + +#endif diff --git a/Modules/QtWidgetsExt/CMakeLists.txt b/Modules/QtWidgetsExt/CMakeLists.txt index 1e0be18163..a8e2d54331 100644 --- a/Modules/QtWidgetsExt/CMakeLists.txt +++ b/Modules/QtWidgetsExt/CMakeLists.txt @@ -1,5 +1,5 @@ MITK_CREATE_MODULE( INCLUDE_DIRS QmitkApplicationBase QmitkPropertyObservers QmitkFunctionalityComponents - DEPENDS MitkAlgorithmsExt MitkQtWidgets + DEPENDS MitkImageStatistics MitkQtWidgets PACKAGE_DEPENDS Qt4|QtWebKit Qwt Qxt ) diff --git a/Modules/QtWidgetsExt/QmitkHistogramJSWidget.cpp b/Modules/QtWidgetsExt/QmitkHistogramJSWidget.cpp index d5dee0a402..cf26a1d79b 100644 --- a/Modules/QtWidgetsExt/QmitkHistogramJSWidget.cpp +++ b/Modules/QtWidgetsExt/QmitkHistogramJSWidget.cpp @@ -1,303 +1,238 @@ /*=================================================================== 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 "QmitkHistogramJSWidget.h" #include "mitkPixelTypeMultiplex.h" #include +#include #include "mitkRenderingManager.h" #include "mitkBaseRenderer.h" #include "mitkImageTimeSelector.h" #include "mitkExtractSliceFilter.h" QmitkHistogramJSWidget::QmitkHistogramJSWidget(QWidget *parent) : QWebView(parent) { // set histogram type to barchart in first instance m_UseLineGraph = false; m_Page = new QmitkJSWebPage(this); setPage(m_Page); // set html from source connect(page()->mainFrame(), SIGNAL(javaScriptWindowObjectCleared()), this, SLOT(AddJSObject())); QUrl myUrl = QUrl("qrc:/qmitk/Histogram.html"); setUrl(myUrl); // set Scrollbars to be always disabled page()->mainFrame()->setScrollBarPolicy(Qt::Horizontal, Qt::ScrollBarAlwaysOff); page()->mainFrame()->setScrollBarPolicy(Qt::Vertical, Qt::ScrollBarAlwaysOff); m_ParametricPath = ParametricPathType::New(); } QmitkHistogramJSWidget::~QmitkHistogramJSWidget() { } // adds an Object of Type QmitkHistogramJSWidget to the JavaScript, using QtWebkitBridge void QmitkHistogramJSWidget::AddJSObject() { page()->mainFrame()->addToJavaScriptWindowObject(QString("histogramData"), this); } // reloads WebView, everytime its size has been changed, so the size of the Histogram fits to the size of the widget void QmitkHistogramJSWidget::resizeEvent(QResizeEvent* resizeEvent) { QWebView::resizeEvent(resizeEvent); // workaround for Qt Bug: https://bugs.webkit.org/show_bug.cgi?id=75984 page()->mainFrame()->evaluateJavaScript("disconnectSignals()"); this->reload(); } // method to expose data to JavaScript by using properties void QmitkHistogramJSWidget::ComputeHistogram(HistogramType* histogram) { m_Histogram = histogram; HistogramConstIteratorType startIt = m_Histogram->End(); HistogramConstIteratorType endIt = m_Histogram->End(); HistogramConstIteratorType it = m_Histogram->Begin(); ClearData(); unsigned int i = 0; bool firstValue = false; // removes frequencies of 0, which are outside the first and last bin for (; it != m_Histogram->End(); ++it) { if (it.GetFrequency() > 0.0) { endIt = it; if (!firstValue) { firstValue = true; startIt = it; } } } ++endIt; // generating Lists of measurement and frequencies for (it = startIt ; it != endIt; ++it, ++i) { QVariant frequency = QVariant::fromValue(it.GetFrequency()); QVariant measurement = it.GetMeasurementVector()[0]; m_Frequency.insert(i, frequency); m_Measurement.insert(i, measurement); } m_IntensityProfile = false; this->SignalDataChanged(); } void QmitkHistogramJSWidget::ClearData() { m_Frequency.clear(); m_Measurement.clear(); } void QmitkHistogramJSWidget::ClearHistogram() { this->ClearData(); this->SignalDataChanged(); } QList QmitkHistogramJSWidget::GetFrequency() { return m_Frequency; } QList QmitkHistogramJSWidget::GetMeasurement() { return m_Measurement; } bool QmitkHistogramJSWidget::GetUseLineGraph() { return m_UseLineGraph; } void QmitkHistogramJSWidget::OnBarRadioButtonSelected() { if (m_UseLineGraph) { m_UseLineGraph = false; this->SignalGraphChanged(); } } void QmitkHistogramJSWidget::OnLineRadioButtonSelected() { if (!m_UseLineGraph) { m_UseLineGraph = true; this->SignalGraphChanged(); } } void QmitkHistogramJSWidget::SetImage(mitk::Image* image) { m_Image = image; } void QmitkHistogramJSWidget::SetPlanarFigure(const mitk::PlanarFigure* planarFigure) { m_PlanarFigure = planarFigure; } template void ReadPixel(mitk::PixelType ptype, mitk::Image::Pointer image, mitk::Index3D indexPoint, double& value) { if (image->GetDimension() == 2) { mitk::ImagePixelReadAccessor readAccess(image, image->GetSliceData(0)); itk::Index<2> idx; idx[0] = indexPoint[0]; idx[1] = indexPoint[1]; value = readAccess.GetPixelByIndex(idx); } else if (image->GetDimension() == 3) { mitk::ImagePixelReadAccessor readAccess(image, image->GetVolumeData(0)); itk::Index<3> idx; idx[0] = indexPoint[0]; idx[1] = indexPoint[1]; idx[2] = indexPoint[2]; value = readAccess.GetPixelByIndex(idx); } else { //unhandled } } void QmitkHistogramJSWidget::ComputeIntensityProfile(unsigned int timeStep) { this->ClearData(); m_ParametricPath->Initialize(); if (m_PlanarFigure.IsNull()) { mitkThrow() << "PlanarFigure not set!"; } if (m_Image.IsNull()) { mitkThrow() << "Image not set!"; } - // Get 2D geometry frame of PlanarFigure - mitk::Geometry2D* planarFigureGeometry2D = dynamic_cast(m_PlanarFigure->GetGeometry(0)); - if (planarFigureGeometry2D == NULL) - { - mitkThrow() << "PlanarFigure has no valid geometry!"; - } - - // Get 3D geometry from Image (needed for conversion of point to index) - mitk::Geometry3D* imageGeometry = m_Image->GetGeometry(0); - if (imageGeometry == NULL) - { - mitkThrow() << "Image has no valid geometry!"; - } - - // Get first poly-line of PlanarFigure (other possible poly-lines in PlanarFigure - // are not supported) - const VertexContainerType vertexContainer = m_PlanarFigure->GetPolyLine(0); + mitk::Image::Pointer image; - VertexContainerType::const_iterator it; - for (it = vertexContainer.begin(); it != vertexContainer.end(); ++it) + if (m_Image->GetDimension() == 4) { - // Map PlanarFigure 2D point to 3D point - mitk::Point3D point3D; - planarFigureGeometry2D->Map(it->Point, point3D); - - // Convert world to index coordinates - mitk::Point3D indexPoint3D; - imageGeometry->WorldToIndex(point3D, indexPoint3D); - - ParametricPathType::OutputType index; - index[0] = indexPoint3D[0]; - index[1] = indexPoint3D[1]; - index[2] = indexPoint3D[2]; + mitk::ImageTimeSelector::Pointer timeSelector = mitk::ImageTimeSelector::New(); + timeSelector->SetInput(m_Image); + timeSelector->SetTimeNr(timeStep); + timeSelector->Update(); - // Add index to parametric path - m_ParametricPath->AddVertex(index); + image = timeSelector->GetOutput(); } - - m_DerivedPath = m_ParametricPath; - - if (m_DerivedPath.IsNull()) + else { - mitkThrow() << "No path set!"; + image = m_Image; } - // Fill item model with line profile data - double distance = 0.0; - mitk::Point3D currentWorldPoint; - - double t; - unsigned int i = 0; - for (i = 0, t = m_DerivedPath->StartOfInput(); ;++i) - { - const PathType::OutputType &continousIndex = m_DerivedPath->Evaluate(t); - - mitk::Point3D worldPoint; - imageGeometry->IndexToWorld(continousIndex, worldPoint); + mitk::IntensityProfile::Pointer intensityProfile = mitk::ComputeIntensityProfile(image, const_cast(m_PlanarFigure.GetPointer())); - if (i == 0) - { - currentWorldPoint = worldPoint; - } - - - distance += currentWorldPoint.EuclideanDistanceTo(worldPoint); - mitk::Index3D indexPoint; - imageGeometry->WorldToIndex(worldPoint, indexPoint); - const mitk::PixelType ptype = m_Image->GetPixelType(); - double intensity = 0.0; - if (m_Image->GetDimension() == 4) - { - mitk::ImageTimeSelector::Pointer timeSelector = mitk::ImageTimeSelector::New(); - timeSelector->SetInput(m_Image); - timeSelector->SetTimeNr(timeStep); - timeSelector->Update(); - mitk::Image::Pointer image = timeSelector->GetOutput(); - mitkPixelTypeMultiplex3( ReadPixel, ptype, image, indexPoint, intensity); - } - else - { - mitkPixelTypeMultiplex3( ReadPixel, ptype, m_Image, indexPoint, intensity); - } + m_Frequency.clear(); + m_Measurement.clear(); - m_Measurement.insert(i, distance); - m_Frequency.insert(i, intensity); + int i = -1; + mitk::IntensityProfile::ConstIterator end = intensityProfile->End(); - // Go to next index; when iteration offset reaches zero, iteration is finished - PathType::OffsetType offset = m_DerivedPath->IncrementInput(t); - if (!(offset[0] || offset[1] || offset[2])) - { - break; - } - - currentWorldPoint = worldPoint; + for (mitk::IntensityProfile::ConstIterator it = intensityProfile->Begin(); it != end; ++it) + { + m_Frequency.push_back(it.GetMeasurementVector()[0]); + m_Measurement.push_back(++i); } m_IntensityProfile = true; m_UseLineGraph = true; this->SignalDataChanged(); } bool QmitkHistogramJSWidget::GetIntensityProfile() { return m_IntensityProfile; }