diff --git a/Modules/Classification/CLCore/src/mitkAbstractGlobalImageFeature.cpp b/Modules/Classification/CLCore/src/mitkAbstractGlobalImageFeature.cpp index 5ac8784d60..ca0d167275 100644 --- a/Modules/Classification/CLCore/src/mitkAbstractGlobalImageFeature.cpp +++ b/Modules/Classification/CLCore/src/mitkAbstractGlobalImageFeature.cpp @@ -1,393 +1,391 @@ /*=================================================================== 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 static void ExtractSlicesFromImages(mitk::Image::Pointer image, mitk::Image::Pointer mask, int direction, std::vector &imageVector, std::vector &maskVector) { typedef itk::Image< double, 2 > FloatImage2DType; typedef itk::Image< unsigned short, 2 > MaskImage2DType; typedef itk::Image< double, 3 > FloatImageType; typedef itk::Image< unsigned short, 3 > MaskImageType; FloatImageType::Pointer itkFloat = FloatImageType::New(); MaskImageType::Pointer itkMask = MaskImageType::New(); mitk::CastToItkImage(mask, itkMask); mitk::CastToItkImage(image, itkFloat); int idxA, idxB, idxC; switch (direction) { case 0: idxA = 1; idxB = 2; idxC = 0; break; case 1: idxA = 0; idxB = 2; idxC = 1; break; case 2: idxA = 0; idxB = 1; idxC = 2; break; default: idxA = 1; idxB = 2; idxC = 0; break; } auto imageSize = image->GetLargestPossibleRegion().GetSize(); FloatImageType::IndexType index3D; FloatImage2DType::IndexType index2D; FloatImage2DType::SpacingType spacing2D; spacing2D[0] = itkFloat->GetSpacing()[idxA]; spacing2D[1] = itkFloat->GetSpacing()[idxB]; for (unsigned int i = 0; i < imageSize[idxC]; ++i) { FloatImage2DType::RegionType region; FloatImage2DType::IndexType start; FloatImage2DType::SizeType size; start[0] = 0; start[1] = 0; size[0] = imageSize[idxA]; size[1] = imageSize[idxB]; region.SetIndex(start); region.SetSize(size); FloatImage2DType::Pointer image2D = FloatImage2DType::New(); image2D->SetRegions(region); image2D->Allocate(); MaskImage2DType::Pointer mask2D = MaskImage2DType::New(); mask2D->SetRegions(region); mask2D->Allocate(); unsigned long voxelsInMask = 0; for (unsigned int a = 0; a < imageSize[idxA]; ++a) { for (unsigned int b = 0; b < imageSize[idxB]; ++b) { index3D[idxA] = a; index3D[idxB] = b; index3D[idxC] = i; index2D[0] = a; index2D[1] = b; image2D->SetPixel(index2D, itkFloat->GetPixel(index3D)); mask2D->SetPixel(index2D, itkMask->GetPixel(index3D)); voxelsInMask += (itkMask->GetPixel(index3D) > 0) ? 1 : 0; } } image2D->SetSpacing(spacing2D); mask2D->SetSpacing(spacing2D); mitk::Image::Pointer tmpFloatImage = mitk::Image::New(); tmpFloatImage->InitializeByItk(image2D.GetPointer()); mitk::GrabItkImageMemory(image2D, tmpFloatImage); mitk::Image::Pointer tmpMaskImage = mitk::Image::New(); tmpMaskImage->InitializeByItk(mask2D.GetPointer()); mitk::GrabItkImageMemory(mask2D, tmpMaskImage); if (voxelsInMask > 0) { imageVector.push_back(tmpFloatImage); maskVector.push_back(tmpMaskImage); } } } std::vector mitk::AbstractGlobalImageFeature::SplitDouble(std::string str, char delimiter) { std::vector internal; std::stringstream ss(str); // Turn the string into a stream. std::string tok; double val; while (std::getline(ss, tok, delimiter)) { std::stringstream s2(tok); s2 >> val; internal.push_back(val); } return internal; } void mitk::AbstractGlobalImageFeature::AddQuantifierArguments(mitkCommandLineParser &parser) { std::string name = GetOptionPrefix(); parser.addArgument(name + "::minimum", name + "::min", mitkCommandLineParser::Float, "Minium Intensity for Quantification", "Defines the minimum Intensity used for Quantification", us::Any()); parser.addArgument(name + "::maximum", name + "::max", mitkCommandLineParser::Float, "Maximum Intensity for Quantification", "Defines the maximum Intensity used for Quantification", us::Any()); parser.addArgument(name + "::bins", name + "::bins", mitkCommandLineParser::Int, "Number of Bins", "Define the number of bins that is used ", us::Any()); parser.addArgument(name + "::binsize", name + "::binsize", mitkCommandLineParser::Float, "Binsize", "Define the size of the used bins", us::Any()); parser.addArgument(name + "::ignore-global-histogram", name + "::ignore-global-histogram", mitkCommandLineParser::Bool, "Ignore the global histogram Parameters", "Ignores the global histogram parameters", us::Any()); parser.addArgument(name + "::ignore-mask-for-histogram", name + "::ignore-mask", mitkCommandLineParser::Bool, "Ignore the global histogram Parameters", "Ignores the global histogram parameters", us::Any()); } void mitk::AbstractGlobalImageFeature::InitializeQuantifierFromParameters(const Image::Pointer & feature, const Image::Pointer &mask, unsigned int defaultBins) { unsigned int bins = 0; double binsize = 0; double minimum = 0; double maximum = 0; auto parsedArgs = GetParameter(); std::string name = GetOptionPrefix(); bool useGlobal = true; if (parsedArgs.count(name + "::ignore-global-histogram")) { useGlobal = false; SetUseMinimumIntensity(false); SetUseMaximumIntensity(false); SetUseBinsize(false); SetUseBins(false); } if (useGlobal) { if (parsedArgs.count("ignore-mask-for-histogram")) { bool tmp = us::any_cast(parsedArgs["ignore-mask-for-histogram"]); SetIgnoreMask(tmp); } if (parsedArgs.count("minimum-intensity")) { minimum = us::any_cast(parsedArgs["minimum-intensity"]); SetMinimumIntensity(minimum); SetUseMinimumIntensity(true); } if (parsedArgs.count("maximum-intensity")) { maximum = us::any_cast(parsedArgs["maximum-intensity"]); SetMaximumIntensity(maximum); SetUseMaximumIntensity(true); } if (parsedArgs.count("bins")) { bins = us::any_cast(parsedArgs["bins"]); SetBins(bins); SetUseBins(true); } if (parsedArgs.count("binsize")) { binsize = us::any_cast(parsedArgs["binsize"]); SetBinsize(binsize); SetUseBinsize(true); } } if (parsedArgs.count(name+"::ignore-mask-for-histogram")) { bool tmp = us::any_cast(parsedArgs[name+"::ignore-mask-for-histogram"]); SetIgnoreMask(tmp); } if (parsedArgs.count(name + "::minimum")) { minimum = us::any_cast(parsedArgs[name + "::minimum"]); SetMinimumIntensity(minimum); SetUseMinimumIntensity(true); } if (parsedArgs.count(name + "::maximum")) { maximum = us::any_cast(parsedArgs[name + "::maximum"]); SetMaximumIntensity(maximum); SetUseMaximumIntensity(true); } if (parsedArgs.count(name + "::bins")) { bins = us::any_cast(parsedArgs[name + "::bins"]); SetBins(bins); } if (parsedArgs.count(name + "::binsize")) { binsize = us::any_cast(parsedArgs[name + "::binsize"]); SetBinsize(binsize); SetUseBinsize(true); } InitializeQuantifier(feature, mask, defaultBins); } void mitk::AbstractGlobalImageFeature::InitializeQuantifier(const Image::Pointer & feature, const Image::Pointer &mask, unsigned int defaultBins) { - MITK_INFO << GetUseMinimumIntensity() << " " << GetUseMaximumIntensity() << " " << GetUseBins() << " " << GetUseBinsize(); - m_Quantifier = IntensityQuantifier::New(); if (GetUseMinimumIntensity() && GetUseMaximumIntensity() && GetUseBinsize()) m_Quantifier->InitializeByBinsizeAndMaximum(GetMinimumIntensity(), GetMaximumIntensity(), GetBinsize()); else if (GetUseMinimumIntensity() && GetUseBins() && GetUseBinsize()) m_Quantifier->InitializeByBinsizeAndBins(GetMinimumIntensity(), GetBins(), GetBinsize()); else if (GetUseMinimumIntensity() && GetUseMaximumIntensity() && GetUseBins()) m_Quantifier->InitializeByMinimumMaximum(GetMinimumIntensity(), GetMaximumIntensity(), GetBins()); // Intialize from Image and Binsize else if (GetUseBinsize() && GetIgnoreMask() && GetUseMinimumIntensity()) m_Quantifier->InitializeByImageAndBinsizeAndMinimum(feature, GetMinimumIntensity(), GetBinsize()); else if (GetUseBinsize() && GetIgnoreMask() && GetUseMaximumIntensity()) m_Quantifier->InitializeByImageAndBinsizeAndMaximum(feature, GetMaximumIntensity(), GetBinsize()); else if (GetUseBinsize() && GetIgnoreMask()) m_Quantifier->InitializeByImageAndBinsize(feature, GetBinsize()); // Initialize form Image, Mask and Binsize else if (GetUseBinsize() && GetUseMinimumIntensity()) m_Quantifier->InitializeByImageRegionAndBinsizeAndMinimum(feature, mask, GetMinimumIntensity(), GetBinsize()); else if (GetUseBinsize() && GetUseMaximumIntensity()) m_Quantifier->InitializeByImageRegionAndBinsizeAndMaximum(feature, mask, GetMaximumIntensity(), GetBinsize()); else if (GetUseBinsize()) m_Quantifier->InitializeByImageRegionAndBinsize(feature, mask, GetBinsize()); // Intialize from Image and Bins else if (GetUseBins() && GetIgnoreMask() && GetUseMinimumIntensity()) m_Quantifier->InitializeByImageAndMinimum(feature, GetMinimumIntensity(), GetBins()); else if (GetUseBins() && GetIgnoreMask() && GetUseMaximumIntensity()) m_Quantifier->InitializeByImageAndMaximum(feature, GetMaximumIntensity(), GetBins()); else if (GetUseBins()) m_Quantifier->InitializeByImage(feature, GetBins()); // Intialize from Image, Mask and Bins else if (GetUseBins() && GetUseMinimumIntensity()) m_Quantifier->InitializeByImageRegionAndMinimum(feature, mask, GetMinimumIntensity(), GetBins()); else if (GetUseBins() && GetUseMaximumIntensity()) m_Quantifier->InitializeByImageRegionAndMaximum(feature, mask, GetMaximumIntensity(), GetBins()); else if (GetUseBins()) m_Quantifier->InitializeByImageRegion(feature, mask, GetBins()); // Default else if (GetIgnoreMask()) m_Quantifier->InitializeByImage(feature, GetBins()); else m_Quantifier->InitializeByImageRegion(feature, mask, defaultBins); } std::string mitk::AbstractGlobalImageFeature::GetCurrentFeatureEncoding() { return ""; } std::string mitk::AbstractGlobalImageFeature::FeatureDescriptionPrefix() { std::string output; output = m_FeatureClassName + "::"; if (m_EncodeParameters) { output += GetCurrentFeatureEncoding() + "::"; } return output; } std::string mitk::AbstractGlobalImageFeature::QuantifierParameterString() { std::stringstream ss; if (GetUseMinimumIntensity() && GetUseMaximumIntensity() && GetUseBinsize()) ss << "Min-" << GetMinimumIntensity() << "_Max-" << GetMaximumIntensity() << "_BS-" << GetBinsize(); else if (GetUseMinimumIntensity() && GetUseBins() && GetUseBinsize()) ss << "Min-" << GetMinimumIntensity() << "_Bins-" << GetBins() << "_BS-" << GetBinsize(); else if (GetUseMinimumIntensity() && GetUseMaximumIntensity() && GetUseBins()) ss << "Min-" << GetMinimumIntensity() << "_Max-" << GetMaximumIntensity() << "_Bins-" << GetBins(); // Intialize from Image and Binsize else if (GetUseBinsize() && GetIgnoreMask() && GetUseMinimumIntensity()) ss << "Min-" << GetMinimumIntensity() << "_BS-" << GetBinsize() << "_FullImage"; else if (GetUseBinsize() && GetIgnoreMask() && GetUseMaximumIntensity()) ss << "Max-" << GetMaximumIntensity() << "_BS-" << GetBinsize() << "_FullImage"; else if (GetUseBinsize() && GetIgnoreMask()) ss << "BS-" << GetBinsize() << "_FullImage"; // Initialize form Image, Mask and Binsize else if (GetUseBinsize() && GetUseMinimumIntensity()) ss << "Min-" << GetMinimumIntensity() << "_BS-" << GetBinsize(); else if (GetUseBinsize() && GetUseMaximumIntensity()) ss << "Max-" << GetMaximumIntensity() << "_BS-" << GetBinsize(); else if (GetUseBinsize()) ss << "BS-" << GetBinsize(); // Intialize from Image and Bins else if (GetUseBins() && GetIgnoreMask() && GetUseMinimumIntensity()) ss << "Min-" << GetMinimumIntensity() << "_Bins-" << GetBins() << "_FullImage"; else if (GetUseBins() && GetIgnoreMask() && GetUseMaximumIntensity()) ss << "Max-" << GetMaximumIntensity() << "_Bins-" << GetBins() << "_FullImage"; else if (GetUseBins()) ss << "Bins-" << GetBins() << "_FullImage"; // Intialize from Image, Mask and Bins else if (GetUseBins() && GetUseMinimumIntensity()) ss << "Min-" << GetMinimumIntensity() << "_Bins-" << GetBins(); else if (GetUseBins() && GetUseMaximumIntensity()) ss << "Max-" << GetMaximumIntensity() << "_Bins-" << GetBins(); else if (GetUseBins()) ss << "Bins-" << GetBins(); // Default else if (GetIgnoreMask()) ss << "Bins-" << GetBins() << "_FullImage"; else ss << "Bins-" << GetBins(); return ss.str(); } mitk::AbstractGlobalImageFeature::FeatureListType mitk::AbstractGlobalImageFeature::CalculateFeaturesSlicewise(const Image::Pointer & feature, const Image::Pointer &mask, int sliceID) { std::vector imageVector; std::vector maskVector; ExtractSlicesFromImages(feature, mask,sliceID, imageVector, maskVector); std::vector statVector; for (std::size_t index = 0; index < imageVector.size(); ++index) { auto stat = this->CalculateFeatures(imageVector[index], maskVector[index]); statVector.push_back(stat); } if (statVector.size() < 1) return FeatureListType(); FeatureListType statMean, statStd, result; for (std::size_t i = 0; i < statVector[0].size(); ++i) { auto cElement1 = statVector[0][i]; cElement1.first = "SliceWise Mean " + cElement1.first; cElement1.second = 0.0; auto cElement2 = statVector[0][i]; cElement2.first = "SliceWise Var. " + cElement2.first; cElement2.second = 0.0; statMean.push_back(cElement1); statStd.push_back(cElement2); } for (auto cStat : statVector) { for (std::size_t i = 0; i < cStat.size(); ++i) { statMean[i].second += cStat[i].second / (1.0*statVector.size()); } } for (auto cStat : statVector) { for (std::size_t i = 0; i < cStat.size(); ++i) { statStd[i].second += (cStat[i].second - statMean[i].second)*(cStat[i].second - statMean[i].second) / (1.0*statVector.size()); } } for (auto cStat : statVector) { std::copy(cStat.begin(), cStat.end(), std::back_inserter(result)); } std::copy(statMean.begin(), statMean.end(), std::back_inserter(result)); std::copy(statStd.begin(), statStd.end(), std::back_inserter(result)); return result; } \ No newline at end of file diff --git a/Modules/Classification/CLUtilities/include/itkLocalIntensity.h b/Modules/Classification/CLUtilities/include/itkLocalIntensity.h new file mode 100644 index 0000000000..186b8a535c --- /dev/null +++ b/Modules/Classification/CLUtilities/include/itkLocalIntensity.h @@ -0,0 +1,176 @@ +/*=================================================================== + +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 itkLocalIntensityFilter_h +#define itkLocalIntensityFilter_h + +#include "itkImageToImageFilter.h" +#include "itkNumericTraits.h" +#include "itkArray.h" +#include "itkSimpleDataObjectDecorator.h" + +namespace itk +{ + + template< typename TInputImage > + class ITK_TEMPLATE_EXPORT LocalIntensityFilter : + public ImageToImageFilter< TInputImage, TInputImage > + { + public: + /** Standard Self typedef */ + typedef LocalIntensityFilter Self; + typedef ImageToImageFilter< TInputImage, TInputImage > Superclass; + typedef SmartPointer< Self > Pointer; + typedef SmartPointer< const Self > ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(LocalIntensityFilter, ImageToImageFilter); + + /** Image related typedefs. */ + typedef typename TInputImage::Pointer InputImagePointer; + + typedef typename TInputImage::RegionType RegionType; + typedef typename TInputImage::SizeType SizeType; + typedef typename TInputImage::IndexType IndexType; + typedef typename TInputImage::PixelType PixelType; + + /** Image related typedefs. */ + itkStaticConstMacro(ImageDimension, unsigned int, + TInputImage::ImageDimension); + + /** Type to use for computations. */ + typedef typename NumericTraits< PixelType >::RealType RealType; + + /** Smart Pointer type to a DataObject. */ + typedef typename DataObject::Pointer DataObjectPointer; + + /** Type of DataObjects used for scalar outputs */ + typedef SimpleDataObjectDecorator< RealType > RealObjectType; + typedef SimpleDataObjectDecorator< PixelType > PixelObjectType; + + /** Return the computed Minimum. */ + PixelType GetMinimum() const + { + return this->GetMinimumOutput()->Get(); + } + PixelObjectType * GetMinimumOutput(); + + const PixelObjectType * GetMinimumOutput() const; + + /** Return the computed Maximum. */ + PixelType GetMaximum() const + { + return this->GetMaximumOutput()->Get(); + } + PixelObjectType * GetMaximumOutput(); + + const PixelObjectType * GetMaximumOutput() const; + + /** Return the computed Mean. */ + RealType GetMean() const + { + return this->GetMeanOutput()->Get(); + } + RealObjectType * GetMeanOutput(); + + const RealObjectType * GetMeanOutput() const; + + /** Return the computed Standard Deviation. */ + RealType GetSigma() const + { + return this->GetSigmaOutput()->Get(); + } + RealObjectType * GetSigmaOutput(); + + const RealObjectType * GetSigmaOutput() const; + + /** Return the computed Variance. */ + RealType GetVariance() const + { + return this->GetVarianceOutput()->Get(); + } + RealObjectType * GetVarianceOutput(); + + const RealObjectType * GetVarianceOutput() const; + + /** Return the compute Sum. */ + RealType GetSum() const + { + return this->GetSumOutput()->Get(); + } + RealObjectType * GetSumOutput(); + + const RealObjectType * GetSumOutput() const; + + /** Make a DataObject of the correct type to be used as the specified + * output. */ + typedef ProcessObject::DataObjectPointerArraySizeType DataObjectPointerArraySizeType; + using Superclass::MakeOutput; + virtual DataObjectPointer MakeOutput(DataObjectPointerArraySizeType idx) ITK_OVERRIDE; + +#ifdef ITK_USE_CONCEPT_CHECKING + // Begin concept checking + itkConceptMacro(InputHasNumericTraitsCheck, + (Concept::HasNumericTraits< PixelType >)); + // End concept checking +#endif + + protected: + LocalIntensityFilter(); + ~LocalIntensityFilter() ITK_OVERRIDE {} + void PrintSelf(std::ostream & os, Indent indent) const ITK_OVERRIDE; + + /** Pass the input through unmodified. Do this by Grafting in the + * AllocateOutputs method. + */ + void AllocateOutputs() ITK_OVERRIDE; + + /** Initialize some accumulators before the threads run. */ + void BeforeThreadedGenerateData() ITK_OVERRIDE; + + /** Do final mean and variance computation from data accumulated in threads. + */ + void AfterThreadedGenerateData() ITK_OVERRIDE; + + /** Multi-thread version GenerateData. */ + void ThreadedGenerateData(const RegionType & + outputRegionForThread, + ThreadIdType threadId) ITK_OVERRIDE; + + // Override since the filter needs all the data for the algorithm + void GenerateInputRequestedRegion() ITK_OVERRIDE; + + // Override since the filter produces all of its output + void EnlargeOutputRequestedRegion(DataObject *data) ITK_OVERRIDE; + + private: + ITK_DISALLOW_COPY_AND_ASSIGN(LocalIntensityFilter); + + Array< RealType > m_ThreadLocalMaximum; + Array< RealType > m_ThreadLocalPeakValue; + Array< RealType > m_ThreadGlobalPeakValue; + }; // end of class +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkLocalIntensityFilter.hxx" +#endif + +#endif + diff --git a/Modules/Classification/CLUtilities/include/itkLocalIntensity.hxx b/Modules/Classification/CLUtilities/include/itkLocalIntensity.hxx new file mode 100644 index 0000000000..94ab5112f3 --- /dev/null +++ b/Modules/Classification/CLUtilities/include/itkLocalIntensity.hxx @@ -0,0 +1,381 @@ +/*=================================================================== + +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 itkLocalIntensityFilter_cpp +#define itkLocalIntensityFilter_cpp + +#include + +#include +#include +#include + +#include + +#include "itkImageScanlineIterator.h" +#include "itkProgressReporter.h" + +namespace itk +{ + template< typename TInputImage > + LocalIntensityFilter< TInputImage > + ::LocalIntensityFilter() :m_ThreadSum(1), m_SumOfSquares(1), m_Count(1), m_ThreadMin(1), m_ThreadMax(1) + { + // first output is a copy of the image, DataObject created by + // superclass + // + // allocate the data objects for the outputs which are + // just decorators around pixel types + for (int i = 1; i < 3; ++i) + { + typename PixelObjectType::Pointer output = + static_cast< PixelObjectType * >(this->MakeOutput(i).GetPointer()); + this->ProcessObject::SetNthOutput(i, output.GetPointer()); + } + // allocate the data objects for the outputs which are + // just decorators around real types + for (int i = 3; i < 7; ++i) + { + typename RealObjectType::Pointer output = + static_cast< RealObjectType * >(this->MakeOutput(i).GetPointer()); + this->ProcessObject::SetNthOutput(i, output.GetPointer()); + } + + this->GetMinimumOutput()->Set(NumericTraits< PixelType >::max()); + this->GetMaximumOutput()->Set(NumericTraits< PixelType >::NonpositiveMin()); + this->GetMeanOutput()->Set(NumericTraits< RealType >::max()); + this->GetSigmaOutput()->Set(NumericTraits< RealType >::max()); + this->GetVarianceOutput()->Set(NumericTraits< RealType >::max()); + this->GetSumOutput()->Set(NumericTraits< RealType >::ZeroValue()); + } + + template< typename TInputImage > + DataObject::Pointer + LocalIntensityFilter< TInputImage > + ::MakeOutput(DataObjectPointerArraySizeType output) + { + switch (output) + { + case 0: + return TInputImage::New().GetPointer(); + break; + case 1: + return PixelObjectType::New().GetPointer(); + break; + case 2: + return PixelObjectType::New().GetPointer(); + break; + case 3: + case 4: + case 5: + case 6: + return RealObjectType::New().GetPointer(); + break; + default: + // might as well make an image + return TInputImage::New().GetPointer(); + break; + } + } + + template< typename TInputImage > + typename LocalIntensityFilter< TInputImage >::PixelObjectType * + LocalIntensityFilter< TInputImage > + ::GetMinimumOutput() + { + return static_cast< PixelObjectType * >(this->ProcessObject::GetOutput(1)); + } + + template< typename TInputImage > + const typename LocalIntensityFilter< TInputImage >::PixelObjectType * + LocalIntensityFilter< TInputImage > + ::GetMinimumOutput() const + { + return static_cast< const PixelObjectType * >(this->ProcessObject::GetOutput(1)); + } + + template< typename TInputImage > + typename LocalIntensityFilter< TInputImage >::PixelObjectType * + LocalIntensityFilter< TInputImage > + ::GetMaximumOutput() + { + return static_cast< PixelObjectType * >(this->ProcessObject::GetOutput(2)); + } + + template< typename TInputImage > + const typename LocalIntensityFilter< TInputImage >::PixelObjectType * + LocalIntensityFilter< TInputImage > + ::GetMaximumOutput() const + { + return static_cast< const PixelObjectType * >(this->ProcessObject::GetOutput(2)); + } + + template< typename TInputImage > + typename LocalIntensityFilter< TInputImage >::RealObjectType * + LocalIntensityFilter< TInputImage > + ::GetMeanOutput() + { + return static_cast< RealObjectType * >(this->ProcessObject::GetOutput(3)); + } + + template< typename TInputImage > + const typename LocalIntensityFilter< TInputImage >::RealObjectType * + LocalIntensityFilter< TInputImage > + ::GetMeanOutput() const + { + return static_cast< const RealObjectType * >(this->ProcessObject::GetOutput(3)); + } + + template< typename TInputImage > + typename LocalIntensityFilter< TInputImage >::RealObjectType * + LocalIntensityFilter< TInputImage > + ::GetSigmaOutput() + { + return static_cast< RealObjectType * >(this->ProcessObject::GetOutput(4)); + } + + template< typename TInputImage > + const typename LocalIntensityFilter< TInputImage >::RealObjectType * + LocalIntensityFilter< TInputImage > + ::GetSigmaOutput() const + { + return static_cast< const RealObjectType * >(this->ProcessObject::GetOutput(4)); + } + + template< typename TInputImage > + typename LocalIntensityFilter< TInputImage >::RealObjectType * + LocalIntensityFilter< TInputImage > + ::GetVarianceOutput() + { + return static_cast< RealObjectType * >(this->ProcessObject::GetOutput(5)); + } + + template< typename TInputImage > + const typename LocalIntensityFilter< TInputImage >::RealObjectType * + LocalIntensityFilter< TInputImage > + ::GetVarianceOutput() const + { + return static_cast< const RealObjectType * >(this->ProcessObject::GetOutput(5)); + } + + template< typename TInputImage > + typename LocalIntensityFilter< TInputImage >::RealObjectType * + LocalIntensityFilter< TInputImage > + ::GetSumOutput() + { + return static_cast< RealObjectType * >(this->ProcessObject::GetOutput(6)); + } + + template< typename TInputImage > + const typename LocalIntensityFilter< TInputImage >::RealObjectType * + LocalIntensityFilter< TInputImage > + ::GetSumOutput() const + { + return static_cast< const RealObjectType * >(this->ProcessObject::GetOutput(6)); + } + + template< typename TInputImage > + void + LocalIntensityFilter< TInputImage > + ::GenerateInputRequestedRegion() + { + Superclass::GenerateInputRequestedRegion(); + if (this->GetInput()) + { + InputImagePointer image = + const_cast< typename Superclass::InputImageType * >(this->GetInput()); + image->SetRequestedRegionToLargestPossibleRegion(); + } + } + + template< typename TInputImage > + void + LocalIntensityFilter< TInputImage > + ::EnlargeOutputRequestedRegion(DataObject *data) + { + Superclass::EnlargeOutputRequestedRegion(data); + data->SetRequestedRegionToLargestPossibleRegion(); + } + + template< typename TInputImage > + void + LocalIntensityFilter< TInputImage > + ::AllocateOutputs() + { + // Pass the input through as the output + InputImagePointer image = + const_cast< TInputImage * >(this->GetInput()); + + this->GraftOutput(image); + + // Nothing that needs to be allocated for the remaining outputs + } + + template< typename TInputImage > + void + LocalIntensityFilter< TInputImage > + ::BeforeThreadedGenerateData() + { + ThreadIdType numberOfThreads = this->GetNumberOfThreads(); + + // Resize the thread temporaries + m_Count.SetSize(numberOfThreads); + m_SumOfSquares.SetSize(numberOfThreads); + m_ThreadSum.SetSize(numberOfThreads); + m_ThreadMin.SetSize(numberOfThreads); + m_ThreadMax.SetSize(numberOfThreads); + + // Initialize the temporaries + m_Count.Fill(NumericTraits< SizeValueType >::ZeroValue()); + m_ThreadSum.Fill(NumericTraits< RealType >::ZeroValue()); + m_SumOfSquares.Fill(NumericTraits< RealType >::ZeroValue()); + m_ThreadMin.Fill(NumericTraits< PixelType >::max()); + m_ThreadMax.Fill(NumericTraits< PixelType >::NonpositiveMin()); + } + + template< typename TInputImage > + void + LocalIntensityFilter< TInputImage > + ::AfterThreadedGenerateData() + { + ThreadIdType i; + SizeValueType count; + RealType sumOfSquares; + + ThreadIdType numberOfThreads = this->GetNumberOfThreads(); + + PixelType minimum; + PixelType maximum; + RealType mean; + RealType sigma; + RealType variance; + RealType sum; + + sum = sumOfSquares = NumericTraits< RealType >::ZeroValue(); + count = 0; + + // Find the min/max over all threads and accumulate count, sum and + // sum of squares + minimum = NumericTraits< PixelType >::max(); + maximum = NumericTraits< PixelType >::NonpositiveMin(); + for (i = 0; i < numberOfThreads; i++) + { + count += m_Count[i]; + sum += m_ThreadSum[i]; + sumOfSquares += m_SumOfSquares[i]; + + if (m_ThreadMin[i] < minimum) + { + minimum = m_ThreadMin[i]; + } + if (m_ThreadMax[i] > maximum) + { + maximum = m_ThreadMax[i]; + } + } + // compute statistics + mean = sum / static_cast< RealType >(count); + + // unbiased estimate + variance = (sumOfSquares - (sum * sum / static_cast< RealType >(count))) + / (static_cast< RealType >(count) - 1); + sigma = std::sqrt(variance); + + // Set the outputs + this->GetMinimumOutput()->Set(minimum); + this->GetMaximumOutput()->Set(maximum); + this->GetMeanOutput()->Set(mean); + this->GetSigmaOutput()->Set(sigma); + this->GetVarianceOutput()->Set(variance); + this->GetSumOutput()->Set(sum); + } + + template< typename TInputImage > + void + LocalIntensityFilter< TInputImage > + ::ThreadedGenerateData(const RegionType & outputRegionForThread, + ThreadIdType threadId) + { + const SizeValueType size0 = outputRegionForThread.GetSize(0); + if (size0 == 0) + { + return; + } + RealType realValue; + PixelType value; + + RealType sum = NumericTraits< RealType >::ZeroValue(); + RealType sumOfSquares = NumericTraits< RealType >::ZeroValue(); + SizeValueType count = NumericTraits< SizeValueType >::ZeroValue(); + PixelType min = NumericTraits< PixelType >::max(); + PixelType max = NumericTraits< PixelType >::NonpositiveMin(); + + ImageScanlineConstIterator< TInputImage > it(this->GetInput(), outputRegionForThread); + + // support progress methods/callbacks + const size_t numberOfLinesToProcess = outputRegionForThread.GetNumberOfPixels() / size0; + ProgressReporter progress(this, threadId, static_cast(numberOfLinesToProcess)); + + // do the work + while (!it.IsAtEnd()) + { + while (!it.IsAtEndOfLine()) + { + value = it.Get(); + realValue = static_cast< RealType >(value); + if (value < min) + { + min = value; + } + if (value > max) + { + max = value; + } + + sum += realValue; + sumOfSquares += (realValue * realValue); + ++count; + ++it; + } + it.NextLine(); + progress.CompletedPixel(); + } + + m_ThreadSum[threadId] = sum; + m_SumOfSquares[threadId] = sumOfSquares; + m_Count[threadId] = count; + m_ThreadMin[threadId] = min; + m_ThreadMax[threadId] = max; + } + + template< typename TImage > + void + LocalIntensityFilter< TImage > + ::PrintSelf(std::ostream & os, Indent indent) const + { + Superclass::PrintSelf(os, indent); + + os << indent << "Minimum: " + << static_cast< typename NumericTraits< PixelType >::PrintType >(this->GetMinimum()) << std::endl; + os << indent << "Maximum: " + << static_cast< typename NumericTraits< PixelType >::PrintType >(this->GetMaximum()) << std::endl; + os << indent << "Sum: " << this->GetSum() << std::endl; + os << indent << "Mean: " << this->GetMean() << std::endl; + os << indent << "Sigma: " << this->GetSigma() << std::endl; + os << indent << "Variance: " << this->GetVariance() << std::endl; + } +} // end namespace itk +#endif \ No newline at end of file diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFLocalIntensity.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFLocalIntensity.cpp index d364e3df15..3e9133a67e 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFLocalIntensity.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFLocalIntensity.cpp @@ -1,179 +1,191 @@ /*=================================================================== 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 // MITK #include #include #include #include #include // ITK #include #include // STL #include struct GIFLocalIntensityParameter { double range; std::string prefix; }; template static void CalculateIntensityPeak(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFLocalIntensity::FeatureListType & featureList, GIFLocalIntensityParameter params) { typedef itk::Image ImageType; typedef itk::Image MaskType; typename MaskType::Pointer itkMask = MaskType::New(); mitk::CastToItkImage(mask, itkMask); double range = params.range; double minimumSpacing = std::numeric_limits::max(); itkImage->GetSpacing(); for (unsigned int i = 0; i < VImageDimension; ++i) { minimumSpacing = (minimumSpacing < itkImage->GetSpacing()[i]) ? minimumSpacing : itkImage->GetSpacing()[i]; } typename ImageType::SizeType regionSize; int offset = std::ceil(range / minimumSpacing); regionSize.Fill(offset); itk::NeighborhoodIterator iter(regionSize, itkImage, itkImage->GetLargestPossibleRegion()); itk::NeighborhoodIterator iterMask(regionSize, itkMask, itkMask->GetLargestPossibleRegion()); typename ImageType::PointType origin; typename ImageType::PointType localPoint; itk::Index index; double tmpPeakValue; double globalPeakValue = 0; double localPeakValue = 0; TPixel localMaximum = 0; + std::vector vectorIsInRange; + index = iter.GetIndex(); + itkImage->TransformIndexToPhysicalPoint(index, origin); + for (itk::SizeValueType i = 0; i < iter.Size(); ++i) + { + itkImage->TransformIndexToPhysicalPoint(iter.GetIndex(i), localPoint); + double dist = origin.EuclideanDistanceTo(localPoint); + vectorIsInRange.push_back((dist < params.range)); + } + int count = 0; + bool inbounds = true; + + iter.NeedToUseBoundaryConditionOff(); + iterMask.NeedToUseBoundaryConditionOff(); + while (!iter.IsAtEnd()) { if (iterMask.GetCenterPixel() > 0) { tmpPeakValue = 0; count = 0; - index = iter.GetIndex(); - itkImage->TransformIndexToPhysicalPoint(index, origin); for (itk::SizeValueType i = 0; i < iter.Size(); ++i) { - itkImage->TransformIndexToPhysicalPoint(iter.GetIndex(i), localPoint); - double dist = origin.EuclideanDistanceTo(localPoint); - if (dist < 6.2) + if (vectorIsInRange[i] ) { if (iter.IndexInBounds(i)) { tmpPeakValue += iter.GetPixel(i); ++count; } } } tmpPeakValue /= count; globalPeakValue = std::max(tmpPeakValue, globalPeakValue); - if (localMaximum == iter.GetCenterPixel()) + auto currentCenterPixelValue = iter.GetCenterPixel(); + if (localMaximum == currentCenterPixelValue) { localPeakValue = std::max(tmpPeakValue,localPeakValue); } - else if (localMaximum < iter.GetCenterPixel()) + else if (localMaximum < currentCenterPixelValue) { - localMaximum = iter.GetCenterPixel(); + localMaximum = currentCenterPixelValue; localPeakValue = tmpPeakValue; } } ++iterMask; ++iter; } featureList.push_back(std::make_pair(params.prefix + "Local Intensity Peak", localPeakValue)); featureList.push_back(std::make_pair(params.prefix + "Global Intensity Peak", globalPeakValue)); } mitk::GIFLocalIntensity::GIFLocalIntensity() : m_Range(6.2) { SetLongName("local-intensity"); SetShortName("loci"); SetFeatureClassName("Local Intensity"); } mitk::GIFLocalIntensity::FeatureListType mitk::GIFLocalIntensity::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; if (image->GetDimension() < 3) { return featureList; } GIFLocalIntensityParameter params; params.range = GetRange(); params.prefix = FeatureDescriptionPrefix(); AccessByItk_3(image, CalculateIntensityPeak, mask, featureList, params); return featureList; } mitk::GIFLocalIntensity::FeatureNameListType mitk::GIFLocalIntensity::GetFeatureNames() { FeatureNameListType featureList; return featureList; } void mitk::GIFLocalIntensity::AddArguments(mitkCommandLineParser &parser) { std::string name = GetOptionPrefix(); parser.addArgument(GetLongName(), name, mitkCommandLineParser::Bool, "Use Local Intensity", "calculates local intensity based features", us::Any()); parser.addArgument(name + "::range", name+"::range", mitkCommandLineParser::Float, "Range for the local intensity", "Give the range that should be used for the local intensity in mm", us::Any()); } void mitk::GIFLocalIntensity::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &, FeatureListType &featureList) { std::string name = GetOptionPrefix(); auto parsedArgs = GetParameter(); if (parsedArgs.count(GetLongName())) { if (parsedArgs.count(name + "::range")) { double range = us::any_cast(parsedArgs[name + "::range"]); this->SetRange(range); } MITK_INFO << "Start calculating local intensity features ...."; auto localResults = this->CalculateFeatures(feature, mask); featureList.insert(featureList.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating local intensity features...."; } } std::string mitk::GIFLocalIntensity::GetCurrentFeatureEncoding() { std::ostringstream ss; ss << m_Range; std::string strRange = ss.str(); return "Range-" + ss.str(); }