diff --git a/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp index 6d33c71c56..9735413843 100644 --- a/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp +++ b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp @@ -1,203 +1,210 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkCommandLineParser.h" #include "mitkImage.h" #include #include #include #include #include #include "mitkIOUtil.h" #include #include #include #include #include "mitkImageAccessByItk.h" #include #include #include #include #include #include #include #include #include #include #include struct statistics_res{ double mean, variance, min, max, count, moment; }; void printstats(statistics_res s) { std::cout << "mean: " << s.mean << std::endl << "variance: " << s.variance << std::endl << "min: " << s.min << std::endl << "max: " << s.max << std::endl << "count: " << s.count << std::endl << "moment: " << s.moment << std::endl; } template void printMap(std::map input) { for (auto it = input.begin(); it != input.end(); ++it) { std::cout << it->first<< ": " << it->second<< std::endl; } std::cout << std::endl; } template < typename TPixel, unsigned int VImageDimension > void get_statistics_boost(itk::Image* itkImage, statistics_res& res){ typedef itk::Image ImageType; itk::ImageRegionConstIterator it(itkImage, itkImage->GetLargestPossibleRegion()); int ctr=0; double sum=0; boost::accumulators::accumulator_set> > acc; for (it.GoToBegin(); !it.IsAtEnd(); ++it) { acc(it.Get()); // currentPixel = it.Get(); // sum+=currentPixel; // ctr+=1; } // res.mean=(double)sum/(double)ctr; res.mean = boost::accumulators::mean(acc); res.variance = boost::accumulators::variance(acc); res.min = boost::accumulators::min(acc); res.max = boost::accumulators::max(acc); res.count = boost::accumulators::count(acc); res.moment = boost::accumulators::moment<2>(acc); std::cout << "sum: " << sum << " N: " << ctr << " mean: " << res.mean << std::endl; } int main( int argc, char* argv[] ) { unsigned int timeStep = 0; std::string inputImageFile; //inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusT.nrrd"; - inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D.nrrd"; + //inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D.nrrd"; + inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/US4DCyl.nrrd"; // Load image mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputImageFile); mitk::Image::Pointer maskImage; std::string maskImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_someSegmentation.nrrd"; maskImage = mitk::IOUtil::LoadImage(maskImageFile); mitk::ImageMaskGenerator::Pointer binaryImageMaskGen = mitk::ImageMaskGenerator::New(); binaryImageMaskGen->SetImageMask(maskImage); mitk::HotspotMaskGenerator::Pointer hotSpotMaskGen = mitk::HotspotMaskGenerator::New(); hotSpotMaskGen->SetImage(inputImage); - hotSpotMaskGen->SetHotspotRadiusInMM(10); + hotSpotMaskGen->SetHotspotRadiusInMM(3); hotSpotMaskGen->SetTimeStep(0); std::string pFfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusT_ellipse.pf"; std::vector loadedObjects; loadedObjects = mitk::IOUtil::Load(pFfile); mitk::BaseData::Pointer pf = loadedObjects[0]; mitk::PlanarFigure::Pointer planarFigure = dynamic_cast(pf.GetPointer()); mitk::PlanarFigureMaskGenerator::Pointer planarFigMaskExtr = mitk::PlanarFigureMaskGenerator::New(); planarFigMaskExtr->SetPlanarFigure(planarFigure); planarFigMaskExtr->SetImage(inputImage); // Calculate statistics mitk::ImageStatisticsCalculator_v2::Pointer calculator = mitk::ImageStatisticsCalculator_v2::New(); std::cout << "calculating statistics (unmasked) itk: " << std::endl; mitk::ImageStatisticsCalculator_v2::StatisticsContainer::Pointer result; calculator->SetInputImage(inputImage); calculator->SetNBinsForHistogramStatistics(100); std::cout << std::endl; // std::cout << "unmasked: " << std::endl; // calculator->SetMask(nullptr); // result = calculator->GetStatistics(timeStep); // result->PrintSelf(); // std::cout << std::endl; // std::cout << "planarfigure: " << std::endl; // calculator->SetMask(planarFigMaskExtr.GetPointer()); // result = calculator->GetStatistics(timeStep, 1); // result->PrintSelf(); // std::cout << std::endl; - std::cout << "ignore pixel value mask: " << std::endl; - mitk::IgnorePixelMaskGenerator::Pointer ignPixVal = mitk::IgnorePixelMaskGenerator::New(); - ignPixVal->SetImage(inputImage); - ignPixVal->SetIgnoredPixelValue(0); - calculator->SetMask(ignPixVal.GetPointer()); - result = calculator->GetStatistics(timeStep, 1); - result->PrintSelf(); - std::cout << std::endl; +// std::cout << "ignore pixel value mask: " << std::endl; +// mitk::IgnorePixelMaskGenerator::Pointer ignPixVal = mitk::IgnorePixelMaskGenerator::New(); +// ignPixVal->SetImage(inputImage); +// ignPixVal->SetIgnoredPixelValue(0); +// calculator->SetMask(ignPixVal.GetPointer()); +// result = calculator->GetStatistics(timeStep, 1); +// result->PrintSelf(); +// std::cout << std::endl; - mitk::IOUtil::SaveImage(ignPixVal->GetMask(), "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_zeroValuedPixelsMask.nrrd"); // std::cout << "image mask: " << std::endl; // calculator->SetMask(binaryImageMaskGen.GetPointer()); // result = calculator->GetStatistics(timeStep, 1); // result->PrintSelf(); -// std::cout << "no mask: " << std::endl; -// calculator->SetMask(nullptr); -// result = calculator->GetStatistics(timeStep); -// result->PrintSelf(); + for (unsigned int i=0; i < inputImage->GetTimeSteps(); i++) + { + std::cout << "no mask: " << std::endl; + std::cout << "timestep: " << i << std::endl; + calculator->SetMask(nullptr); + result = calculator->GetStatistics(i); + result->PrintSelf(); + std::cout << std::endl; + } +// std::cout << "hotspot mask: " << std::endl; // calculator->SetMask(hotSpotMaskGen.GetPointer()); // result = calculator->GetStatistics(timeStep, 1); // result->PrintSelf(); +// mitk::IOUtil::SaveImage(hotSpotMaskGen->GetMask(), "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusTHotspot.nrrd"); // std::cout << std::endl << std::endl << "calculating statistics old imgstatcalc: " << std::endl; // mitk::ImageStatisticsCalculator::Statistics statisticsStruct; // mitk::ImageStatisticsCalculator::Pointer calculator_old = mitk::ImageStatisticsCalculator::New(); // calculator_old->SetImage(inputImage); // calculator_old->SetImageMask(maskImage); // //calculator_old->SetMaskingModeToImage(); // calculator_old->SetHistogramBinSize(100); // calculator_old->SetPlanarFigure(planarFigure); // calculator_old->SetMaskingModeToPlanarFigure(); // calculator_old->SetMaskingModeToNone(); // calculator_old->SetDoIgnorePixelValue(false); // calculator_old->ComputeStatistics(timeStep); // statisticsStruct = calculator_old->GetStatistics(timeStep); // std::cout << "calculating statistics boost: " << std::endl; // statistics_res res; // AccessByItk_n(inputImage, get_statistics_boost, (res)); // printstats(res); return EXIT_SUCCESS; } diff --git a/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp b/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp index bc14376d56..a95eb142fa 100644 --- a/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp @@ -1,535 +1,542 @@ #include #include #include #include #include #include "mitkImageAccessByItk.h" #include #include #include namespace mitk { HotspotMaskGenerator::HotspotMaskGenerator(): m_HotspotRadiusinMM(6.2035049089940), // radius of a 1cm3 sphere in mm m_HotspotMustBeCompletelyInsideImage(true), m_Label(1) { - m_Modified=true; + m_TimeStep = 0; + m_Modified = true; } void HotspotMaskGenerator::SetImage(mitk::Image::Pointer inputImage) { if (inputImage != m_Image) { m_Image = inputImage; m_Modified = true; } } void HotspotMaskGenerator::SetMask(MaskGenerator::Pointer mask) { if (mask != m_Mask) { m_Mask = mask; m_Modified = true; } } HotspotMaskGenerator::~HotspotMaskGenerator() { } void HotspotMaskGenerator::SetHotspotRadiusInMM(double radiusInMillimeter) { if(radiusInMillimeter != m_HotspotRadiusinMM) { m_HotspotRadiusinMM = radiusInMillimeter; m_Modified = true; } } const double& HotspotMaskGenerator::GetHotspotRadiusinMM() const { return m_HotspotRadiusinMM; } bool HotspotMaskGenerator::GetHotspotMustBeCompletelyInsideImage() const { return m_HotspotMustBeCompletelyInsideImage; } mitk::Image::Pointer HotspotMaskGenerator::GetMask() { if (m_Modified) { if ( m_Image.IsNull() ) { throw std::runtime_error( "Error: image empty!" ); } if ( m_TimeStep >= m_Image->GetTimeSteps() ) { throw std::runtime_error( "Error: invalid time step!" ); } mitk::ImageTimeSelector::Pointer imageTimeSelector = mitk::ImageTimeSelector::New(); imageTimeSelector->SetInput( m_Image ); imageTimeSelector->SetTimeNr( m_TimeStep ); imageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image::Pointer timeSliceImage = imageTimeSelector->GetOutput(); m_internalImage = timeSliceImage; m_internalMask2D = nullptr; // is this correct when this variable holds a smart pointer? m_internalMask3D = nullptr; if ( m_Mask != nullptr ) { m_Mask->SetTimeStep(m_TimeStep); mitk::Image::Pointer timeSliceMask = m_Mask->GetMask(); if ( m_internalImage->GetDimension() == 3 ) { CastToItkImage(timeSliceMask, m_internalMask3D); AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D, m_Label); } else if ( m_internalImage->GetDimension() == 2 ) { CastToItkImage(timeSliceMask, m_internalMask2D); AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D, m_Label); } else { throw std::runtime_error( "Error: invalid image dimension" ); } } else { if ( m_internalImage->GetDimension() == 3 ) { AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D, m_Label); } else if ( m_internalImage->GetDimension() == 2 ) { AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D, m_Label); } else { throw std::runtime_error( "Error: invalid image dimension" ); } } m_Modified = false; } return m_InternalMask; } void HotspotMaskGenerator::SetLabel(unsigned int label) { if (label != m_Label) { m_Label = label; m_Modified = true; } } template HotspotMaskGenerator::ImageExtrema HotspotMaskGenerator::CalculateExtremaWorld( const itk::Image* inputImage, typename itk::Image::Pointer maskImage, double neccessaryDistanceToImageBorderInMM, unsigned int label ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef itk::ImageRegionConstIteratorWithIndex MaskImageIteratorType; typedef itk::ImageRegionConstIteratorWithIndex InputImageIndexIteratorType; typename ImageType::SpacingType spacing = inputImage->GetSpacing(); ImageExtrema minMax; minMax.Defined = false; minMax.MaxIndex.set_size(VImageDimension); minMax.MaxIndex.set_size(VImageDimension); typename ImageType::RegionType allowedExtremaRegion = inputImage->GetLargestPossibleRegion(); bool keepDistanceToImageBorders( neccessaryDistanceToImageBorderInMM > 0 ); if (keepDistanceToImageBorders) { long distanceInPixels[VImageDimension]; for(unsigned short dimension = 0; dimension < VImageDimension; ++dimension) { // To confirm that the whole hotspot is inside the image we have to keep a specific distance to the image-borders, which is as long as // the radius. To get the amount of indices we divide the radius by spacing and add 0.5 because voxels are center based: // For example with a radius of 2.2 and a spacing of 1 two indices are enough because 2.2 / 1 + 0.5 = 2.7 => 2. // But with a radius of 2.7 we need 3 indices because 2.7 / 1 + 0.5 = 3.2 => 3 distanceInPixels[dimension] = int( neccessaryDistanceToImageBorderInMM / spacing[dimension] + 0.5); } allowedExtremaRegion.ShrinkByRadius(distanceInPixels); } InputImageIndexIteratorType imageIndexIt(inputImage, allowedExtremaRegion); float maxValue = itk::NumericTraits::min(); float minValue = itk::NumericTraits::max(); typename ImageType::IndexType maxIndex; typename ImageType::IndexType minIndex; for(unsigned short i = 0; i < VImageDimension; ++i) { maxIndex[i] = 0; minIndex[i] = 0; } if (maskImage != nullptr) { MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); typename ImageType::IndexType imageIndex; typename ImageType::PointType worldPosition; typename ImageType::IndexType maskIndex; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { imageIndex = maskIndex = maskIt.GetIndex(); if(maskIt.Get() == label) { if( allowedExtremaRegion.IsInside(imageIndex) ) { imageIndexIt.SetIndex( imageIndex ); double value = imageIndexIt.Get(); minMax.Defined = true; //Calculate minimum, maximum and corresponding index-values if( value > maxValue ) { maxIndex = imageIndexIt.GetIndex(); maxValue = value; } if(value < minValue ) { minIndex = imageIndexIt.GetIndex(); minValue = value; } } } } } else { for(imageIndexIt.GoToBegin(); !imageIndexIt.IsAtEnd(); ++imageIndexIt) { double value = imageIndexIt.Get(); minMax.Defined = true; //Calculate minimum, maximum and corresponding index-values if( value > maxValue ) { maxIndex = imageIndexIt.GetIndex(); maxValue = value; } if(value < minValue ) { minIndex = imageIndexIt.GetIndex(); minValue = value; } } } minMax.MaxIndex.set_size(VImageDimension); minMax.MinIndex.set_size(VImageDimension); for(unsigned int i = 0; i < minMax.MaxIndex.size(); ++i) { minMax.MaxIndex[i] = maxIndex[i]; } for(unsigned int i = 0; i < minMax.MinIndex.size(); ++i) { minMax.MinIndex[i] = minIndex[i]; } minMax.Max = maxValue; minMax.Min = minValue; return minMax; } template itk::Size HotspotMaskGenerator::CalculateConvolutionKernelSize( double spacing[VImageDimension], double radiusInMM ) { typedef itk::Image< float, VImageDimension > KernelImageType; typedef typename KernelImageType::SizeType SizeType; SizeType maskSize; for(unsigned int i = 0; i < VImageDimension; ++i) { maskSize[i] = static_cast( 2 * radiusInMM / spacing[i]); // We always want an uneven size to have a clear center point in the convolution mask if(maskSize[i] % 2 == 0 ) { ++maskSize[i]; } } return maskSize; } template itk::SmartPointer< itk::Image > HotspotMaskGenerator::GenerateHotspotSearchConvolutionKernel(double mmPerPixel[VImageDimension], double radiusInMM ) { std::stringstream ss; for (unsigned int i = 0; i < VImageDimension; ++i) { ss << mmPerPixel[i]; if (i < VImageDimension -1) ss << ","; } MITK_DEBUG << "Update convolution kernel for spacing (" << ss.str() << ") and radius " << radiusInMM << "mm"; double radiusInMMSquared = radiusInMM * radiusInMM; typedef itk::Image< float, VImageDimension > KernelImageType; typename KernelImageType::Pointer convolutionKernel = KernelImageType::New(); // Calculate size and allocate mask image typedef typename KernelImageType::SizeType SizeType; SizeType maskSize = this->CalculateConvolutionKernelSize(mmPerPixel, radiusInMM); mitk::Point3D convolutionMaskCenterIndex; convolutionMaskCenterIndex.Fill(0.0); for(unsigned int i = 0; i < VImageDimension; ++i) { convolutionMaskCenterIndex[i] = 0.5 * (double)(maskSize[i]-1); } typedef typename KernelImageType::IndexType IndexType; IndexType maskIndex; maskIndex.Fill(0); typedef typename KernelImageType::RegionType RegionType; RegionType maskRegion; maskRegion.SetSize(maskSize); maskRegion.SetIndex(maskIndex); convolutionKernel->SetRegions(maskRegion); convolutionKernel->SetSpacing(mmPerPixel); convolutionKernel->Allocate(); // Fill mask image values by subsampling the image grid typedef itk::ImageRegionIteratorWithIndex MaskIteratorType; MaskIteratorType maskIt(convolutionKernel,maskRegion); int numberOfSubVoxelsPerDimension = 2; // per dimension! int numberOfSubVoxels = ::pow( static_cast(numberOfSubVoxelsPerDimension), static_cast(VImageDimension) ); double subVoxelSizeInPixels = 1.0 / (double)numberOfSubVoxelsPerDimension; double valueOfOneSubVoxel = 1.0 / (double)numberOfSubVoxels; double maskValue = 0.0; mitk::Point3D subVoxelIndexPosition; double distanceSquared = 0.0; typedef itk::ContinuousIndex ContinuousIndexType; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { ContinuousIndexType indexPoint(maskIt.GetIndex()); mitk::Point3D voxelPosition; for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) { voxelPosition[dimension] = indexPoint[dimension]; } maskValue = 0.0; mitk::Vector3D subVoxelOffset; subVoxelOffset.Fill(0.0); // iterate sub-voxels by iterating all possible offsets for (subVoxelOffset[0] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[0] < +0.5; subVoxelOffset[0] += subVoxelSizeInPixels) { for (subVoxelOffset[1] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[1] < +0.5; subVoxelOffset[1] += subVoxelSizeInPixels) { for (subVoxelOffset[2] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[2] < +0.5; subVoxelOffset[2] += subVoxelSizeInPixels) { subVoxelIndexPosition = voxelPosition + subVoxelOffset; // this COULD be integrated into the for-loops if neccessary (add voxelPosition to initializer and end condition) distanceSquared = (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] * (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] + (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] * (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] + (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2] * (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2]; if (distanceSquared <= radiusInMMSquared) { maskValue += valueOfOneSubVoxel; } } } } maskIt.Set( maskValue ); } return convolutionKernel; } template itk::SmartPointer > HotspotMaskGenerator::GenerateConvolutionImage( const itk::Image* inputImage ) { double mmPerPixel[VImageDimension]; for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) { mmPerPixel[dimension] = inputImage->GetSpacing()[dimension]; } // update convolution kernel typedef itk::Image< float, VImageDimension > KernelImageType; typename KernelImageType::Pointer convolutionKernel = this->GenerateHotspotSearchConvolutionKernel(mmPerPixel, m_HotspotRadiusinMM); // update convolution image typedef itk::Image< TPixel, VImageDimension > InputImageType; typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; typedef itk::FFTConvolutionImageFilter ConvolutionFilterType; typename ConvolutionFilterType::Pointer convolutionFilter = ConvolutionFilterType::New(); typedef itk::ConstantBoundaryCondition BoundaryConditionType; BoundaryConditionType boundaryCondition; boundaryCondition.SetConstant(0.0); if (m_HotspotMustBeCompletelyInsideImage) { // overwrite default boundary condition convolutionFilter->SetBoundaryCondition(&boundaryCondition); } convolutionFilter->SetInput(inputImage); convolutionFilter->SetKernelImage(convolutionKernel); convolutionFilter->SetNormalize(true); MITK_DEBUG << "Update Convolution image for hotspot search"; convolutionFilter->UpdateLargestPossibleRegion(); typename ConvolutionImageType::Pointer convolutionImage = convolutionFilter->GetOutput(); convolutionImage->SetSpacing( inputImage->GetSpacing() ); // only workaround because convolution filter seems to ignore spacing of input image return convolutionImage; } template < typename TPixel, unsigned int VImageDimension> void HotspotMaskGenerator::FillHotspotMaskPixels( itk::Image* maskImage, itk::Point sphereCenter, double sphereRadiusInMM ) { typedef itk::Image< TPixel, VImageDimension > MaskImageType; typedef itk::ImageRegionIteratorWithIndex MaskImageIteratorType; MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); typename MaskImageType::IndexType maskIndex; typename MaskImageType::PointType worldPosition; // this is not very smart. I would rather use a 0 initialized mask (not the case here -> blame CalculateHotspotMask) and find the region where I need to iterate over, then iterate only over the small region for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { maskIndex = maskIt.GetIndex(); maskImage->TransformIndexToPhysicalPoint(maskIndex, worldPosition); maskIt.Set( worldPosition.EuclideanDistanceTo(sphereCenter) <= sphereRadiusInMM ? 1 : 0 ); } } template void HotspotMaskGenerator::CalculateHotspotMask(itk::Image* inputImage, typename itk::Image::Pointer maskImage, unsigned int label) { typedef itk::Image< TPixel, VImageDimension > InputImageType; typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; typedef itk::Image< float, VImageDimension > KernelImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; - typedef itk::ImageDuplicator< InputImageType > DuplicatorType; - typedef itk::CastImageFilter< InputImageType, MaskImageType > CastFilterType; typename ConvolutionImageType::Pointer convolutionImage = this->GenerateConvolutionImage(inputImage); if (convolutionImage.IsNull()) { MITK_ERROR << "Empty convolution image in CalculateHotspotStatistics(). We should never reach this state (logic error)."; throw std::logic_error("Empty convolution image in CalculateHotspotStatistics()"); } // if mask image is not defined, create an image of the same size as inputImage and fill it with 1's // there is maybe a better way to do this!? if (maskImage == nullptr) { maskImage = MaskImageType::New(); typename MaskImageType::RegionType maskRegion = inputImage->GetLargestPossibleRegion(); typename MaskImageType::SpacingType maskSpacing = inputImage->GetSpacing(); typename MaskImageType::PointType maskOrigin = inputImage->GetOrigin(); typename MaskImageType::DirectionType maskDirection = inputImage->GetDirection(); maskImage->SetRegions(maskRegion); maskImage->Allocate(); maskImage->SetOrigin(maskOrigin); maskImage->SetSpacing(maskSpacing); maskImage->SetDirection(maskDirection); maskImage->FillBuffer(1); label = 1; } // find maximum in convolution image, given the current mask double requiredDistanceToBorder = m_HotspotMustBeCompletelyInsideImage ? m_HotspotRadiusinMM : -1.0; ImageExtrema convolutionImageInformation = CalculateExtremaWorld(convolutionImage.GetPointer(), maskImage, requiredDistanceToBorder, label); bool isHotspotDefined = convolutionImageInformation.Defined; if (!isHotspotDefined) { MITK_ERROR << "No origin of hotspot-sphere was calculated!"; m_InternalMask = nullptr; } else { // create a binary mask around the "hotspot" region, fill the shape of a sphere around our hotspot center - typename DuplicatorType::Pointer copyMachine = DuplicatorType::New(); - copyMachine->SetInputImage(inputImage); - copyMachine->Update(); - - typename CastFilterType::Pointer caster = CastFilterType::New(); - caster->SetInput( copyMachine->GetOutput() ); - caster->Update(); - typename MaskImageType::Pointer hotspotMaskITK = caster->GetOutput(); +// typename DuplicatorType::Pointer copyMachine = DuplicatorType::New(); +// copyMachine->SetInputImage(inputImage); +// copyMachine->Update(); + +// typename CastFilterType::Pointer caster = CastFilterType::New(); +// caster->SetInput( copyMachine->GetOutput() ); +// caster->Update(); + typename MaskImageType::Pointer hotspotMaskITK = MaskImageType::New(); + hotspotMaskITK->SetOrigin(inputImage->GetOrigin()); + hotspotMaskITK->SetSpacing(inputImage->GetSpacing()); + hotspotMaskITK->SetLargestPossibleRegion(inputImage->GetLargestPossibleRegion()); + hotspotMaskITK->SetBufferedRegion(inputImage->GetBufferedRegion()); + hotspotMaskITK->SetDirection(inputImage->GetDirection()); + hotspotMaskITK->SetNumberOfComponentsPerPixel(inputImage->GetNumberOfComponentsPerPixel()); + hotspotMaskITK->Allocate(); + hotspotMaskITK->FillBuffer(1); typedef typename InputImageType::IndexType IndexType; IndexType maskCenterIndex; for (unsigned int d =0; d< VImageDimension;++d) { maskCenterIndex[d]=convolutionImageInformation.MaxIndex[d]; } typename ConvolutionImageType::PointType maskCenter; inputImage->TransformIndexToPhysicalPoint(maskCenterIndex,maskCenter); FillHotspotMaskPixels(hotspotMaskITK.GetPointer(), maskCenter, m_HotspotRadiusinMM); //obtain mitk::Image::Pointer from itk::Image mitk::Image::Pointer hotspotMaskAsMITKImage = mitk::GrabItkImageMemory(hotspotMaskITK); m_InternalMask = hotspotMaskAsMITKImage; } } } diff --git a/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h b/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h index e060ebb244..aa5fff82d5 100644 --- a/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h +++ b/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h @@ -1,53 +1,56 @@ #ifndef MITKIGNOREPIXELMASKGEN_ #define MITKIGNOREPIXELMASKGEN_ #include #include #include #include #include namespace mitk { class MITKIMAGESTATISTICS_EXPORT IgnorePixelMaskGenerator: public MaskGenerator { public: /** Standard Self typedef */ typedef IgnorePixelMaskGenerator Self; typedef MaskGenerator Superclass; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; typedef double RealType; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(IgnorePixelMaskGenerator, MaskGenerator) void SetImage(mitk::Image::Pointer image); void SetIgnoredPixelValue(RealType pixelValue); mitk::Image::Pointer GetMask(); protected: IgnorePixelMaskGenerator(): m_IgnoredPixelValue(std::numeric_limits::min()) - {} + { + m_TimeStep = 0; + m_Modified = true; + } ~IgnorePixelMaskGenerator(){} template void InternalCalculateMask(typename itk::Image* image); private: mitk::Image::Pointer m_Image; RealType m_IgnoredPixelValue; }; } #endif diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp index 99c69df1ab..17fdd8bf26 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp @@ -1,525 +1,516 @@ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "itkImageFileWriter.h" namespace mitk { void ImageStatisticsCalculator_v2::SetInputImage(mitk::Image::Pointer image) { if (image != m_Image) { m_Image = image; m_StatisticsUpdateRequiredByTimeStep.resize(m_Image->GetTimeSteps()); m_StatisticsByTimeStep.resize(m_Image->GetTimeSteps()); this->SetAllStatisticsToUpdateRequired(); } } void ImageStatisticsCalculator_v2::SetMask(mitk::MaskGenerator::Pointer mask) { if (mask != m_MaskGenerator) { m_MaskGenerator = mask; this->SetAllStatisticsToUpdateRequired(); } } void ImageStatisticsCalculator_v2::SetNBinsForHistogramStatistics(unsigned int nBins) { if (nBins != m_nBinsForHistogramStatistics) { m_nBinsForHistogramStatistics = nBins; this->SetAllStatisticsToUpdateRequired(); } } unsigned int ImageStatisticsCalculator_v2::GetNBinsForHistogramStatistics() const { return m_nBinsForHistogramStatistics; } ImageStatisticsCalculator_v2::StatisticsContainer::Pointer ImageStatisticsCalculator_v2::GetStatistics(unsigned int timeStep, unsigned int label) { std::cout << "getting statistics..." << std::endl; if (timeStep >= m_StatisticsByTimeStep.size()) { throw std::runtime_error("invalid timeStep in ImageStatisticsCalculator_v2::GetStatistics"); } if (m_Image.IsNull()) { throw std::runtime_error("no image"); } if (m_MaskGenerator.IsNotNull()) { m_MaskGenerator->SetTimeStep(timeStep); m_InternalMask = m_MaskGenerator->GetMask(); } ImageTimeSelector::Pointer imgTimeSel = ImageTimeSelector::New(); imgTimeSel->SetInput(m_Image); imgTimeSel->SetTimeNr(timeStep); imgTimeSel->UpdateLargestPossibleRegion(); m_ImageTimeSlice = imgTimeSel->GetOutput(); if (m_StatisticsUpdateRequiredByTimeStep[timeStep]) { std::cout << "calculating statistics..." << std::endl; // Calculate statistics with/without mask if (m_MaskGenerator.IsNull()) { // 1) calculate statistics unmasked: // plug image into itkstatisticsimagefilter (will be replaced by my awesome filter later on) // retrieve statistics and save them std::cout << "unmasked" << std::endl; AccessByItk_1(m_ImageTimeSlice, InternalCalculateStatisticsUnmasked, timeStep) } else { std::cout << "masked" << std::endl; // 2) calculate statistics masked // extract mask image region // plug mask and image into itklabelstatisticsimagefilter AccessByItk_1(m_ImageTimeSlice, InternalCalculateStatisticsMasked, timeStep) } m_StatisticsUpdateRequiredByTimeStep[timeStep]=false; } for (std::vector::iterator it = m_StatisticsByTimeStep[timeStep].begin(); it != m_StatisticsByTimeStep[timeStep].end(); ++it) { StatisticsContainer::Pointer statCont = *it; if (statCont->GetLabel() == label) { return statCont; } } MITK_ERROR << "Invalid Label: " << label; } void ImageStatisticsCalculator_v2::SetAllStatisticsToUpdateRequired() { for (unsigned int i = 0; i < m_StatisticsUpdateRequiredByTimeStep.size(); i++) { this->SetStatsTimeStepToUpdateRequired(i); } } void ImageStatisticsCalculator_v2::SetStatsTimeStepToUpdateRequired(unsigned int timeStep) { if (timeStep >= m_StatisticsUpdateRequiredByTimeStep.size()) { throw std::runtime_error("invalid timeStep in ImageStatisticsCalculator_v2::SetStatsTimeStepToUpdateRequired"); } m_StatisticsUpdateRequiredByTimeStep[timeStep] = true; m_StatisticsByTimeStep[timeStep].resize(0); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::InternalCalculateStatisticsUnmasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep) { typedef typename itk::Image< TPixel, VImageDimension > ImageType; StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New(); typename itk::ExtendedStatisticsImageFilter::Pointer statisticsFilter = itk::ExtendedStatisticsImageFilter::New(); statisticsFilter->SetInput(image); statisticsFilter->SetCoordinateTolerance(0.001); statisticsFilter->SetDirectionTolerance(0.001); // TODO: this is single threaded. Implement our own image filter that does this multi threaded typename itk::MinimumMaximumImageCalculator::Pointer imgMinMaxFilter = itk::MinimumMaximumImageCalculator::New(); imgMinMaxFilter->SetImage(image); imgMinMaxFilter->Compute(); vnl_vector minIndex, maxIndex; typename ImageType::IndexType tmpMinIndex = imgMinMaxFilter->GetIndexOfMinimum(); typename ImageType::IndexType tmpMaxIndex = imgMinMaxFilter->GetIndexOfMaximum(); minIndex.set_size(tmpMaxIndex.GetIndexDimension()); maxIndex.set_size(tmpMaxIndex.GetIndexDimension()); for (unsigned int i=0; i < tmpMaxIndex.GetIndexDimension(); i++) { minIndex[i] = tmpMinIndex[i]; maxIndex[i] = tmpMaxIndex[i]; } statisticsResult->SetMinIndex(minIndex); statisticsResult->SetMaxIndex(maxIndex); // TODO: REMOVE THIS, the histogram statistics should be a separate module statisticsFilter->SetBinSize(m_nBinsForHistogramStatistics); try { statisticsFilter->Update(); } catch (const itk::ExceptionObject& e) { mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } // no mask, therefore just one label = the whole image m_StatisticsByTimeStep[timeStep].resize(1); statisticsResult->SetLabel(1); statisticsResult->SetN(image->GetLargestPossibleRegion().GetNumberOfPixels()); statisticsResult->SetMean(statisticsFilter->GetMean()); statisticsResult->SetMin(statisticsFilter->GetMinimum()); statisticsResult->SetMax(statisticsFilter->GetMaximum()); statisticsResult->SetVariance(statisticsFilter->GetVariance()); statisticsResult->SetStd(statisticsFilter->GetSigma()); statisticsResult->SetSkewness(statisticsFilter->GetSkewness()); statisticsResult->SetKurtosis(statisticsFilter->GetKurtosis()); statisticsResult->SetRMS(std::sqrt(std::pow(statisticsFilter->GetMean(), 2.) + statisticsFilter->GetVariance())); // variance = sigma^2 statisticsResult->SetMPP(statisticsFilter->GetMPP()); // histogram statistics HistogramType::Pointer histogram = InternalCalculateHistogramUnmasked(image, statisticsResult->GetMin(), statisticsResult->GetMax()); statisticsResult->SetHistogram(histogram); HistogramStatisticsCalculator histStatCalc; histStatCalc.SetHistogram(histogram); histStatCalc.CalculateStatistics(); statisticsResult->SetEntropy(histStatCalc.GetEntropy()); statisticsResult->SetMedian(histStatCalc.GetMedian()); statisticsResult->SetUniformity(histStatCalc.GetUniformity()); statisticsResult->SetUPP(histStatCalc.GetUPP()); m_StatisticsByTimeStep[timeStep][0] = statisticsResult; } template < typename TPixel, unsigned int VImageDimension > ImageStatisticsCalculator_v2::HistogramType::Pointer ImageStatisticsCalculator_v2::InternalCalculateHistogramUnmasked( typename itk::Image< TPixel, VImageDimension >* image, double minVal, double maxVal) { typedef typename itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Statistics::ImageToHistogramFilter HistogramFilterType; typedef itk::Statistics::Histogram HistogramType; typedef typename HistogramFilterType::HistogramSizeType SizeType; typename HistogramFilterType::Pointer histogramFilter = HistogramFilterType::New(); histogramFilter->SetInput(image); SizeType histSize(1); histSize[0] = m_nBinsForHistogramStatistics; histogramFilter->SetHistogramSize(histSize); // is that correct?? typename HistogramFilterType::HistogramMeasurementVectorType lowerBound(1); typename HistogramFilterType::HistogramMeasurementVectorType upperBound(1); lowerBound[0] = minVal; upperBound[0] = maxVal; histogramFilter->SetHistogramBinMinimum(lowerBound); histogramFilter->SetHistogramBinMaximum(upperBound); histogramFilter->Update(); typename HistogramType::Pointer histogram = histogramFilter->GetOutput(); return histogram; } - /* Taken from original ImageStatisticsCalculator - * This function is not ideal on so many levels - 1) it does not support more than one label - 2) it assumes that the label id is 1 (which does not have to be the case) - 3) single threaded - => I suggest writing an itk filter that finds min and max value. This filter could then be called from our new - statisticsimagefilter in the beforethreadedgenereatedata (provided that a min/max for the histogram is not yet defined) - EDIT: I changed this function so that it simply searches the image min and max and completely disregards the mask. - Since we call this function after we cropped the image to the mask region we can run it on the LargestPossibleRegion of the image */ - template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::GetMinAndMaxValue( - typename itk::Image< TPixel, VImageDimension >* inputImage, - double &minVal, - double &maxVal - ) - { - typedef itk::Image< TPixel, VImageDimension > ImageType; - - typedef itk::ImageRegionConstIteratorWithIndex ImageRegionIteratorType; - - ImageRegionIteratorType imageIterator( inputImage, inputImage->GetLargestPossibleRegion() ); - - maxVal = std::numeric_limits::min(); - minVal = std::numeric_limits::max(); - long counter(0); - double actualPixelValue(0); - - for( imageIterator.GoToBegin(); !imageIterator.IsAtEnd(); ++imageIterator) - { - counter++; - - actualPixelValue = imageIterator.Value(); - - if(actualPixelValue > maxVal) - { - maxVal = actualPixelValue; - } - else if(actualPixelValue < minVal) - { - minVal = actualPixelValue; - } - } - } + /* the whole min and max calculation thing is completely bad... need a new itk filter for that */ +// template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::GetMinAndMaxValue( +// typename itk::Image< TPixel, VImageDimension >* inputImage, +// double &minVal, +// double &maxVal +// ) +// { +// typedef itk::Image< TPixel, VImageDimension > ImageType; + +// typedef itk::ImageRegionConstIteratorWithIndex ImageRegionIteratorType; + +// ImageRegionIteratorType imageIterator( inputImage, inputImage->GetLargestPossibleRegion() ); + +// maxVal = std::numeric_limits::min(); +// minVal = std::numeric_limits::max(); +// long counter(0); +// double actualPixelValue(0); + +// for( imageIterator.GoToBegin(); !imageIterator.IsAtEnd(); ++imageIterator) +// { +// counter++; + +// actualPixelValue = imageIterator.Value(); + +// if(actualPixelValue > maxVal) +// { +// maxVal = actualPixelValue; +// } +// else if(actualPixelValue < minVal) +// { +// minVal = actualPixelValue; +// } +// } +// } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::InternalCalculateStatisticsMasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskType; typedef itk::ExtendedLabelStatisticsImageFilter< ImageType, MaskType > ImageStatisticsFilterType; typedef MaskUtilities< TPixel, VImageDimension > MaskUtilType; typedef itk::MaskImageFilter2 MaskImageFilterType2; typedef itk::MaskImageFilter MaskImageFilterType; typedef itk::MinimumMaximumImageCalculator MinMaxCalculatorType; // maskImage has to have the same dimension as image typename MaskType::Pointer maskImage = MaskType::New(); maskImage = ImageToItkImage< unsigned short, VImageDimension >(m_InternalMask); typename MaskUtilType::Pointer maskUtil = MaskUtilType::New(); maskUtil->SetImage(image); maskUtil->SetMask(maskImage.GetPointer()); // if mask is smaller than image, extract the image region where the mask is typename ImageType::Pointer adaptedImage = ImageType::New(); adaptedImage = maskUtil->ExtractMaskImageRegion(); // this also checks mask sanity // double minVal, maxVal; // GetMinAndMaxValue(adaptedImage.GetPointer(), minVal, maxVal); /** Searching the min and max value is only done over the pixels that are nonzero in the mask. This bevahior differs from the 'old' * imagestatisticscalculator where min and max were looked for over the entire adaptedImage*/ // find min, max, minindex and maxindex // make sure to only look in the masked region, use a masker for this typename MaskImageFilterType::Pointer maskImageFilter = MaskImageFilterType::New(); maskImageFilter->SetInput1(adaptedImage); maskImageFilter->SetInput2(maskImage); maskImageFilter->SetCoordinateTolerance(0.001); maskImageFilter->SetDirectionTolerance(0.001); maskImageFilter->Update(); typename MinMaxCalculatorType::Pointer imgMinMaxFilter = MinMaxCalculatorType::New(); imgMinMaxFilter->SetImage(maskImageFilter->GetOutput()); imgMinMaxFilter->Compute(); double minVal, maxVal; minVal = imgMinMaxFilter->GetMinimum(); maxVal = imgMinMaxFilter->GetMaximum(); typename ImageStatisticsFilterType::Pointer imageStatisticsFilter = ImageStatisticsFilterType::New(); imageStatisticsFilter->SetDirectionTolerance(0.001); imageStatisticsFilter->SetCoordinateTolerance(0.001); imageStatisticsFilter->SetInput(adaptedImage); imageStatisticsFilter->SetLabelInput(maskImage); imageStatisticsFilter->SetHistogramParameters(m_nBinsForHistogramStatistics, floor(minVal), ceil(maxVal)); imageStatisticsFilter->Update(); std::list labels = imageStatisticsFilter->GetRelevantLabels(); std::list::iterator it = labels.begin(); m_StatisticsByTimeStep[timeStep].resize(0); while(it != labels.end()) { StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New(); // find min, max, minindex and maxindex // make sure to only look in the masked region, use a masker for this typename MaskImageFilterType2::Pointer maskImageFilter = MaskImageFilterType2::New(); maskImageFilter->SetInput1(adaptedImage); maskImageFilter->SetInput2(maskImage); maskImageFilter->SetOutsideValue(0); maskImageFilter->SetMaskingValue(*it); - double tmp_label = *it; maskImageFilter->SetCoordinateTolerance(0.001); maskImageFilter->SetDirectionTolerance(0.001); maskImageFilter->UpdateLargestPossibleRegion(); // std::string outfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_masked_image.nrrd"; // typedef itk::ImageFileWriter< ImageType > WriterType; // typename WriterType::Pointer writer = WriterType::New(); // writer->SetFileName(outfile); // writer->SetInput(maskImageFilter->GetOutput()); // writer->Update(); // TODO adapt minindex and maxindex if mask was 2D (careful: a mask can be 2d even if the itk::Image is 3D (3d image with a size of 1 in one dimension)) typename MinMaxCalculatorType::Pointer imgMinMaxFilter = MinMaxCalculatorType::New(); imgMinMaxFilter->SetImage(maskImageFilter->GetOutput()); imgMinMaxFilter->Compute(); vnl_vector minIndex, maxIndex; typename ImageType::IndexType tmpMinIndex = imgMinMaxFilter->GetIndexOfMinimum(); typename ImageType::IndexType tmpMaxIndex = imgMinMaxFilter->GetIndexOfMaximum(); minIndex.set_size(tmpMaxIndex.GetIndexDimension()); maxIndex.set_size(tmpMaxIndex.GetIndexDimension()); for (unsigned int i=0; i < tmpMaxIndex.GetIndexDimension(); i++) { minIndex[i] = tmpMinIndex[i] + (maskImage->GetOrigin()[i] - image->GetOrigin()[i]) / (double) maskImage->GetSpacing()[i]; maxIndex[i] = tmpMaxIndex[i] + (maskImage->GetOrigin()[i] - image->GetOrigin()[i]) / (double) maskImage->GetSpacing()[i]; } statisticsResult->SetMinIndex(minIndex); statisticsResult->SetMaxIndex(maxIndex); statisticsResult->SetN(imageStatisticsFilter->GetSum(*it) / (double) imageStatisticsFilter->GetMean(*it)); statisticsResult->SetMean(imageStatisticsFilter->GetMean(*it)); statisticsResult->SetMin(imageStatisticsFilter->GetMinimum(*it)); statisticsResult->SetMax(imageStatisticsFilter->GetMaximum(*it)); statisticsResult->SetVariance(imageStatisticsFilter->GetVariance(*it)); statisticsResult->SetStd(imageStatisticsFilter->GetSigma(*it)); statisticsResult->SetSkewness(imageStatisticsFilter->GetSkewness(*it)); statisticsResult->SetKurtosis(imageStatisticsFilter->GetKurtosis(*it)); statisticsResult->SetRMS(std::sqrt(std::pow(imageStatisticsFilter->GetMean(*it), 2.) + imageStatisticsFilter->GetVariance(*it))); // variance = sigma^2 statisticsResult->SetMPP(imageStatisticsFilter->GetMPP(*it)); statisticsResult->SetLabel(*it); HistogramType::Pointer histogram = imageStatisticsFilter->GetHistogram(*it); statisticsResult->SetHistogram(histogram); HistogramStatisticsCalculator histStatCalc; histStatCalc.SetHistogram(histogram); histStatCalc.CalculateStatistics(); statisticsResult->SetEntropy(histStatCalc.GetEntropy()); statisticsResult->SetMedian(histStatCalc.GetMedian()); statisticsResult->SetUniformity(histStatCalc.GetUniformity()); statisticsResult->SetUPP(histStatCalc.GetUPP()); m_StatisticsByTimeStep[timeStep].push_back(statisticsResult); ++it; } } ImageStatisticsCalculator_v2::StatisticsContainer::StatisticsContainer(): m_N(std::numeric_limits::max()), m_Mean(std::numeric_limits::max()), m_Min(std::numeric_limits::min()), m_Max(std::numeric_limits::max()), m_Std(std::numeric_limits::max()), m_Variance(std::numeric_limits::max()), m_Skewness(std::numeric_limits::max()), m_Kurtosis(std::numeric_limits::max()), m_RMS(std::numeric_limits::max()), m_MPP(std::numeric_limits::max()), m_Median(std::numeric_limits::max()), m_Uniformity(std::numeric_limits::max()), m_UPP(std::numeric_limits::max()), m_Entropy(std::numeric_limits::max()) { m_minIndex.set_size(0); m_maxIndex.set_size(0); } ImageStatisticsCalculator_v2::statisticsMapType ImageStatisticsCalculator_v2::StatisticsContainer::GetStatisticsAsMap() { ImageStatisticsCalculator_v2::statisticsMapType statisticsAsMap; statisticsAsMap["N"] = m_N; statisticsAsMap["Mean"] = m_Mean; statisticsAsMap["Min"] = m_Min; statisticsAsMap["Max"] = m_Max; statisticsAsMap["StandardDeviation"] = m_Std; statisticsAsMap["Variance"] = m_Variance; statisticsAsMap["Skewness"] = m_Skewness; statisticsAsMap["Kurtosis"] = m_Kurtosis; statisticsAsMap["RMS"] = m_RMS; statisticsAsMap["MPP"] = m_MPP; statisticsAsMap["Median"] = m_Median; statisticsAsMap["Uniformity"] = m_Uniformity; statisticsAsMap["UPP"] = m_UPP; statisticsAsMap["Entropy"] = m_Entropy; return statisticsAsMap; } void ImageStatisticsCalculator_v2::StatisticsContainer::Reset() { m_N = std::numeric_limits::max(); m_Mean = std::numeric_limits::max(); m_Min = std::numeric_limits::max(); m_Max = std::numeric_limits::max(); m_Std = std::numeric_limits::max(); m_Variance = std::numeric_limits::max(); m_Skewness = std::numeric_limits::max(); m_Kurtosis = std::numeric_limits::max(); m_RMS = std::numeric_limits::max(); m_MPP = std::numeric_limits::max(); m_Median = std::numeric_limits::max(); m_Uniformity = std::numeric_limits::max(); m_UPP = std::numeric_limits::max(); m_Entropy = std::numeric_limits::max(); } void ImageStatisticsCalculator_v2::StatisticsContainer::PrintSelf() { ImageStatisticsCalculator_v2::statisticsMapType statMap = this->GetStatisticsAsMap(); // print all map key value pairs for (auto it = statMap.begin(); it != statMap.end(); ++it) { std::cout << it->first << ": " << it->second << std::endl; } // print the min and max index std::cout << "Min Index:" << std::endl; for (auto it = this->GetMinIndex().begin(); it != this->GetMinIndex().end(); it++) { std::cout << *it << " "; } std::cout << std::endl; // print the min and max index std::cout << "Max Index:" << std::endl; for (auto it = this->GetMaxIndex().begin(); it != this->GetMaxIndex().end(); it++) { std::cout << *it << " "; } std::cout << std::endl; } } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h index 9eecd3d6aa..fdd8af3207 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h @@ -1,333 +1,333 @@ #ifndef MITKIMAGESTATISTICSCALCULATOR2 #define MITKIMAGESTATISTICSCALCULATOR2 #include #include #include #include #include #include #include /*#define mitkGetRefConstMacro(name, type) \ const type& Get##name() const \ { \ return &m_##name; \ } \ \ type& Get##name() \ { \ return &m_##name; \ } \ */ namespace mitk { class MITKIMAGESTATISTICS_EXPORT ImageStatisticsCalculator_v2: public itk::Object { public: /** Standard Self typedef */ typedef ImageStatisticsCalculator_v2 Self; typedef itk::Object Superclass; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(ImageStatisticsCalculator_v2, itk::Object) typedef double statisticsValueType; typedef std::map statisticsMapType; typedef itk::Statistics::Histogram HistogramType; class StatisticsContainer: public itk::Object { public: /** Standard Self typedef */ typedef StatisticsContainer Self; typedef itk::Object Superclass; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(StatisticsContainer, itk::Object) typedef double RealType; statisticsMapType GetStatisticsAsMap(); void Reset(); void SetN(long n) { m_N = n; } const long& GetN() const { return m_N; } void SetMean(RealType mean) { m_Mean = mean; } const RealType& GetMean() const { return m_Mean; } void SetVariance(RealType variance) { m_Variance = variance; } const RealType& GetVariance() const { return m_Variance; } void SetStd(RealType std) { m_Std = std; } const RealType& GetStd() const { return m_Std; } void SetMin(RealType minVal) { m_Min = minVal; } const RealType& GetMin() const { return m_Min; } void SetMax(RealType maxVal) { m_Max = maxVal; } const RealType& GetMax() const { return m_Max; } void SetRMS(RealType rms) { m_RMS = rms; } const RealType& GetRMS() const { return m_RMS; } void SetSkewness(RealType skewness) { m_Skewness = skewness; } const RealType& GetSkewness() const { return m_Skewness; } void SetKurtosis(RealType kurtosis) { m_Kurtosis = kurtosis; } const RealType& GetKurtosis() const { return m_Kurtosis; } void SetMPP(RealType mpp) { m_MPP = mpp; } const RealType& GetMPP() const { return m_MPP; } void SetLabel(unsigned int label) { m_Label = label; } const unsigned int& GetLabel() const { return m_Label; } void SetMinIndex(vnl_vector& minIndex) { m_minIndex = minIndex; } const vnl_vector& GetMinIndex() const { return m_minIndex; } void SetMaxIndex(vnl_vector& maxIndex) { m_maxIndex = maxIndex; } const vnl_vector& GetMaxIndex() const { return m_maxIndex; } void SetHistogram(HistogramType::Pointer hist) { if (m_Histogram != hist) { m_Histogram = hist; } } const HistogramType::Pointer GetHistogram() const { return m_Histogram; } void SetEntropy(RealType entropy) { m_Entropy = entropy; } const RealType & GetEntropy() const { return m_Entropy; } void SetMedian(RealType median) { m_Median = median; } const RealType & GetMedian() const { return m_Median; } void SetUniformity(RealType uniformity) { m_Uniformity = uniformity; } const RealType & GetUniformity() const { return m_Uniformity; } void SetUPP(RealType upp) { m_UPP = upp; } const RealType & GetUPP() const { return m_UPP; } void PrintSelf(); protected: StatisticsContainer(); private: // not pretty, is temporary long m_N; RealType m_Mean, m_Min, m_Max, m_Std, m_Variance; RealType m_Skewness; RealType m_Kurtosis; RealType m_RMS; RealType m_MPP; vnl_vector m_minIndex, m_maxIndex; RealType m_Median; RealType m_Uniformity; RealType m_UPP; RealType m_Entropy; unsigned int m_Label; HistogramType::Pointer m_Histogram; }; void SetInputImage(mitk::Image::Pointer image); void SetMask(mitk::MaskGenerator::Pointer mask); void SetNBinsForHistogramStatistics(unsigned int nBins); unsigned int GetNBinsForHistogramStatistics() const; StatisticsContainer::Pointer GetStatistics(unsigned int timeStep=0, unsigned int label=1); protected: ImageStatisticsCalculator_v2(){ m_nBinsForHistogramStatistics = 100; }; private: void SetAllStatisticsToUpdateRequired(); void SetStatsTimeStepToUpdateRequired(unsigned int timeStep); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsUnmasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep); template < typename TPixel, unsigned int VImageDimension > typename HistogramType::Pointer InternalCalculateHistogramUnmasked( typename itk::Image< TPixel, VImageDimension >* image, double minVal, double maxVal); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsMasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep); - template < typename TPixel, unsigned int VImageDimension > void GetMinAndMaxValue( - typename itk::Image< TPixel, VImageDimension >* inputImage, - double &minVal, - double &maxVal - ); +// template < typename TPixel, unsigned int VImageDimension > void GetMinAndMaxValue( +// typename itk::Image< TPixel, VImageDimension >* inputImage, +// double &minVal, +// double &maxVal +// ); std::string GetNameOfClass() { return std::string("ImageStatisticsCalculator_v2"); } mitk::Image::Pointer m_Image; mitk::Image::Pointer m_ImageTimeSlice; mitk::MaskGenerator::Pointer m_MaskGenerator; mitk::Image::Pointer m_InternalMask; unsigned int m_nBinsForHistogramStatistics; std::vector m_StatisticsUpdateRequiredByTimeStep; // holds which time steps are valid and which ones have to be recalculated std::vector> m_StatisticsByTimeStep; }; } #endif // MITKIMAGESTATISTICSCALCULATOR2 diff --git a/Modules/ImageStatistics/mitkMaskGenerator.h b/Modules/ImageStatistics/mitkMaskGenerator.h index 0a863852fd..74dca4a577 100644 --- a/Modules/ImageStatistics/mitkMaskGenerator.h +++ b/Modules/ImageStatistics/mitkMaskGenerator.h @@ -1,48 +1,60 @@ #ifndef MITKMASKGENERATOR #define MITKMASKGENERATOR #include #include #include #include #include namespace mitk { +/** +* \class MaskGenerator +* \brief Base Class for all Mask Generators. Mask generators are classes that provide functionality for the +* creation of binary (or unsigned short) masks that can be applied to an image. See dervied classes for more +* information. +*/ class MITKIMAGESTATISTICS_EXPORT MaskGenerator: public itk::Object { public: /** Standard Self typedef */ typedef MaskGenerator Self; typedef itk::Object Superclass; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(MaskGenerator, itk::Object) //~MaskGenerator(); + /** + * @brief GetMask must be overridden by derived classes. + * @return mitk::Image::Pointer of generated mask + */ virtual mitk::Image::Pointer GetMask(); + /** + * @brief SetTimeStep is used to set the time step for which the mask is to be generated + * @param timeStep + */ void SetTimeStep(unsigned int timeStep); protected: MaskGenerator(); - void CheckMaskSanity(mitk::Image::Pointer image); - unsigned int m_TimeStep; mitk::Image::Pointer m_InternalMask; bool m_Modified; private: }; } #endif // MITKMASKGENERATOR diff --git a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp index 2c2554d8bd..4266f175b2 100644 --- a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp @@ -1,402 +1,404 @@ #include #include #include #include "mitkImageAccessByItk.h" #include #include #include #include #include #include #include #include #include #include #include #include #include namespace mitk { void PlanarFigureMaskGenerator::SetPlanarFigure(mitk::PlanarFigure::Pointer planarFigure) { if ( planarFigure.IsNull() ) { throw std::runtime_error( "Error: planar figure empty!" ); } if ( !planarFigure->IsClosed() ) { throw std::runtime_error( "Masking not possible for non-closed figures" ); } const PlaneGeometry *planarFigurePlaneGeometry = planarFigure->GetPlaneGeometry(); if ( planarFigurePlaneGeometry == nullptr ) { throw std::runtime_error( "Planar-Figure not yet initialized!" ); } const PlaneGeometry *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry ); if ( planarFigureGeometry == nullptr ) { throw std::runtime_error( "Non-planar planar figures not supported!" ); } if (planarFigure != m_PlanarFigure) { m_Modified = true; m_PlanarFigure = planarFigure; } } void PlanarFigureMaskGenerator::SetImage(mitk::Image::Pointer image) { // check dimension unsigned int dimension = image->GetDimension(); if (dimension > 3) { MITK_ERROR << "unsupported image dimension"; } const BaseGeometry *imageGeometry = image->GetGeometry(); if ( imageGeometry == nullptr ) { throw std::runtime_error( "Image geometry invalid!" ); } if (image != m_InternalInputImage) { m_Modified = true; m_InternalInputImage = image; } } template < typename TPixel, unsigned int VImageDimension > void PlanarFigureMaskGenerator::InternalCalculateMaskFromPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, 2 > MaskImage2DType; - typedef itk::CastImageFilter< ImageType, MaskImage2DType > CastFilterType; +// typedef itk::CastImageFilter< ImageType, MaskImage2DType > CastFilterType; // Generate mask image as new image with same header as input image and // initialize with 1. - typename CastFilterType::Pointer castFilter = CastFilterType::New(); - castFilter->SetInput( image ); - castFilter->Update(); - castFilter->GetOutput()->FillBuffer( 1 ); -// typename itk::Image::Pointer maskImage = itk::Image::New(); -// maskImage->SetDirection(image->GetDirection()); -// maskImage->SetOrigin(image->GetOrigin()); -// maskImage->SetSpacing(image->GetSpacing()); -// maskImage->SetLargestPossibleRegion(image->GetLargestPossibleRegion()); -// maskImage->Allocate(); -// maskImage->FillBuffer(1); +// typename CastFilterType::Pointer castFilter = CastFilterType::New(); +// castFilter->SetInput( image ); +// castFilter->Update(); +// castFilter->GetOutput()->FillBuffer( 1 ); + typename MaskImage2DType::Pointer maskImage = MaskImage2DType::New(); + maskImage->SetOrigin(image->GetOrigin()); + maskImage->SetSpacing(image->GetSpacing()); + maskImage->SetLargestPossibleRegion(image->GetLargestPossibleRegion()); + maskImage->SetBufferedRegion(image->GetBufferedRegion()); + maskImage->SetDirection(image->GetDirection()); + maskImage->SetNumberOfComponentsPerPixel(image->GetNumberOfComponentsPerPixel()); + maskImage->Allocate(); + maskImage->FillBuffer(1); // all PolylinePoints of the PlanarFigure are stored in a vtkPoints object. // These points are used by the vtkLassoStencilSource to create // a vtkImageStencil. const mitk::PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); const mitk::BaseGeometry *imageGeometry3D = m_InternalInputImage->GetGeometry( 0 ); // If there is a second poly line in a closed planar figure, treat it as a hole. PlanarFigure::PolyLineType planarFigureHolePolyline; if (m_PlanarFigure->GetPolyLinesSize() == 2) planarFigureHolePolyline = m_PlanarFigure->GetPolyLine(1); // Determine x- and y-dimensions depending on principal axis // TODO use plane geometry normal to determine that automatically, then check whether the PF is aligned with one of the three principal axis int i0, i1; switch ( axis ) { case 0: i0 = 1; i1 = 2; break; case 1: i0 = 0; i1 = 2; break; case 2: default: i0 = 0; i1 = 1; break; } // store the polyline contour as vtkPoints object bool outOfBounds = false; vtkSmartPointer points = vtkSmartPointer::New(); typename PlanarFigure::PolyLineType::const_iterator it; for ( it = planarFigurePolyline.begin(); it != planarFigurePolyline.end(); ++it ) { Point3D point3D; // Convert 2D point back to the local index coordinates of the selected // image // Fabian: From PlaneGeometry documentation: // Converts a 2D point given in mm (pt2d_mm) relative to the upper-left corner of the geometry into the corresponding world-coordinate (a 3D point in mm, pt3d_mm). // To convert a 2D point given in units (e.g., pixels in case of an image) into a 2D point given in mm (as required by this method), use IndexToWorld. planarFigurePlaneGeometry->Map( *it, point3D ); // Polygons (partially) outside of the image bounds can not be processed // further due to a bug in vtkPolyDataToImageStencil if ( !imageGeometry3D->IsInside( point3D ) ) { outOfBounds = true; } imageGeometry3D->WorldToIndex( point3D, point3D ); points->InsertNextPoint( point3D[i0], point3D[i1], 0 ); } vtkSmartPointer holePoints = nullptr; if (!planarFigureHolePolyline.empty()) { holePoints = vtkSmartPointer::New(); Point3D point3D; PlanarFigure::PolyLineType::const_iterator end = planarFigureHolePolyline.end(); for (it = planarFigureHolePolyline.begin(); it != end; ++it) { // Fabian: same as above planarFigurePlaneGeometry->Map(*it, point3D); imageGeometry3D->WorldToIndex(point3D, point3D); holePoints->InsertNextPoint(point3D[i0], point3D[i1], 0); } } // mark a malformed 2D planar figure ( i.e. area = 0 ) as out of bounds // this can happen when all control points of a rectangle lie on the same line = two of the three extents are zero double bounds[6] = {0, 0, 0, 0, 0, 0}; points->GetBounds( bounds ); bool extent_x = (fabs(bounds[0] - bounds[1])) < mitk::eps; bool extent_y = (fabs(bounds[2] - bounds[3])) < mitk::eps; bool extent_z = (fabs(bounds[4] - bounds[5])) < mitk::eps; // throw an exception if a closed planar figure is deformed, i.e. has only one non-zero extent if ( m_PlanarFigure->IsClosed() && ((extent_x && extent_y) || (extent_x && extent_z) || (extent_y && extent_z))) { mitkThrow() << "Figure has a zero area and cannot be used for masking."; } if ( outOfBounds ) { throw std::runtime_error( "Figure at least partially outside of image bounds!" ); } // create a vtkLassoStencilSource and set the points of the Polygon vtkSmartPointer lassoStencil = vtkSmartPointer::New(); lassoStencil->SetShapeToPolygon(); lassoStencil->SetPoints( points ); vtkSmartPointer holeLassoStencil = nullptr; if (holePoints.GetPointer() != nullptr) { holeLassoStencil = vtkSmartPointer::New(); holeLassoStencil->SetShapeToPolygon(); holeLassoStencil->SetPoints(holePoints); } // Export from ITK to VTK (to use a VTK filter) typedef itk::VTKImageImport< MaskImage2DType > ImageImportType; typedef itk::VTKImageExport< MaskImage2DType > ImageExportType; typename ImageExportType::Pointer itkExporter = ImageExportType::New(); - //itkExporter->SetInput( maskImage ); - itkExporter->SetInput( castFilter->GetOutput() ); + itkExporter->SetInput( maskImage ); +// itkExporter->SetInput( castFilter->GetOutput() ); vtkSmartPointer vtkImporter = vtkSmartPointer::New(); this->ConnectPipelines( itkExporter, vtkImporter ); // Apply the generated image stencil to the input image vtkSmartPointer imageStencilFilter = vtkSmartPointer::New(); imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); imageStencilFilter->SetStencilConnection(lassoStencil->GetOutputPort()); imageStencilFilter->ReverseStencilOff(); imageStencilFilter->SetBackgroundValue( 0 ); imageStencilFilter->Update(); vtkSmartPointer holeStencilFilter = nullptr; if (holeLassoStencil.GetPointer() != nullptr) { holeStencilFilter = vtkSmartPointer::New(); holeStencilFilter->SetInputConnection(imageStencilFilter->GetOutputPort()); holeStencilFilter->SetStencilConnection(holeLassoStencil->GetOutputPort()); holeStencilFilter->ReverseStencilOn(); holeStencilFilter->SetBackgroundValue(0); holeStencilFilter->Update(); } // Export from VTK back to ITK vtkSmartPointer vtkExporter = vtkSmartPointer::New(); vtkExporter->SetInputConnection( holeStencilFilter.GetPointer() == nullptr ? imageStencilFilter->GetOutputPort() : holeStencilFilter->GetOutputPort()); vtkExporter->Update(); typename ImageImportType::Pointer itkImporter = ImageImportType::New(); this->ConnectPipelines( vtkExporter, itkImporter ); itkImporter->Update(); typedef itk::ImageDuplicator< ImageImportType::OutputImageType > DuplicatorType; DuplicatorType::Pointer duplicator = DuplicatorType::New(); duplicator->SetInputImage( itkImporter->GetOutput() ); duplicator->Update(); // Store mask m_InternalITKImageMask2D = duplicator->GetOutput(); } bool PlanarFigureMaskGenerator::GetPrincipalAxis( const BaseGeometry *geometry, Vector3D vector, unsigned int &axis ) { vector.Normalize(); for ( unsigned int i = 0; i < 3; ++i ) { Vector3D axisVector = geometry->GetAxisVector( i ); axisVector.Normalize(); if ( fabs( fabs( axisVector * vector ) - 1.0) < mitk::eps ) { axis = i; return true; } } return false; } void PlanarFigureMaskGenerator::CalculateMask() { if (m_InternalInputImage->GetTimeSteps() > 0) { mitk::ImageTimeSelector::Pointer imgTimeSel = mitk::ImageTimeSelector::New(); imgTimeSel->SetInput(m_InternalInputImage); imgTimeSel->SetTimeNr(m_TimeStep); imgTimeSel->UpdateLargestPossibleRegion(); m_InternalTimeSliceImage = imgTimeSel->GetOutput(); } else { m_InternalTimeSliceImage = m_InternalInputImage; } m_InternalITKImageMask2D = nullptr; const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); const PlaneGeometry *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry ); const BaseGeometry *imageGeometry = m_InternalInputImage->GetGeometry(); // Find principal direction of PlanarFigure in input image unsigned int axis; if ( !this->GetPrincipalAxis( imageGeometry, planarFigureGeometry->GetNormal(), axis ) ) { throw std::runtime_error( "Non-aligned planar figures not supported!" ); } m_PlanarFigureAxis = axis; // Find slice number corresponding to PlanarFigure in input image typename itk::Image< unsigned short, 3 >::IndexType index; imageGeometry->WorldToIndex( planarFigureGeometry->GetOrigin(), index ); unsigned int slice = index[axis]; // extract image slice which corresponds to the planarFigure and sotre it in m_InternalImageSlice mitk::Image::Pointer inputImageSlice = extract2DImageSlice(axis, slice); // Compute mask from PlanarFigure AccessFixedDimensionByItk_1(inputImageSlice, InternalCalculateMaskFromPlanarFigure, 2, axis) //convert itk mask to mitk::Image::Pointer and return it mitk::Image::Pointer planarFigureMaskImage; planarFigureMaskImage = mitk::GrabItkImageMemory(m_InternalITKImageMask2D); planarFigureMaskImage->SetGeometry(inputImageSlice->GetGeometry()); Convert2Dto3DImageFilter::Pointer sliceTo3DImageConverter = Convert2Dto3DImageFilter::New(); sliceTo3DImageConverter->SetInput(planarFigureMaskImage); sliceTo3DImageConverter->Update(); m_InternalMask = sliceTo3DImageConverter->GetOutput(); } mitk::Image::Pointer PlanarFigureMaskGenerator::GetMask() { if (m_Modified) { if (m_InternalInputImage.IsNull()) { MITK_ERROR << "Image is not set."; } if (m_PlanarFigure.IsNull()) { MITK_ERROR << "PlanarFigure is not set."; } if (m_TimeStep != 0) { MITK_WARN << "Multiple TimeSteps are not supported in PlanarFigureMaskGenerator (yet)."; } this->CalculateMask(); } m_Modified = false; return m_InternalMask; } mitk::Image::Pointer PlanarFigureMaskGenerator::extract2DImageSlice(unsigned int axis, unsigned int slice) { // Extract slice with given position and direction from image unsigned int dimension = m_InternalTimeSliceImage->GetDimension(); mitk::Image::Pointer imageSlice = mitk::Image::New(); if (dimension == 3) { ExtractImageFilter::Pointer imageExtractor = ExtractImageFilter::New(); imageExtractor->SetInput( m_InternalTimeSliceImage ); imageExtractor->SetSliceDimension( axis ); imageExtractor->SetSliceIndex( slice ); imageExtractor->Update(); imageSlice = imageExtractor->GetOutput(); } else if(dimension == 2) { imageSlice = m_InternalTimeSliceImage; } else { MITK_ERROR << "Unsupported image dimension. Dimension is: " << dimension << ". Only 2D and 3D images are supported."; } return imageSlice; } } diff --git a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h index 48e56b7aa1..e8bcd4b177 100644 --- a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h +++ b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h @@ -1,107 +1,118 @@ #ifndef MITKPLANARFIGUREMASKGENERATOR #define MITKPLANARFIGUREMASKGENERATOR #include #include #include #include #include #include #include #include #include #include namespace mitk { +/** +* \class PlanarFigureMaskGenerator +* \brief Derived from MaskGenerator. This class is used to convert a mitk::PlanarFigure into a binary image mask +*/ class MITKIMAGESTATISTICS_EXPORT PlanarFigureMaskGenerator: public MaskGenerator { public: /** Standard Self typedef */ typedef PlanarFigureMaskGenerator Self; typedef MaskGenerator Superclass; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(PlanarFigureMaskGenerator, MaskGenerator) - mitk::Image::Pointer GetMask(); + /** + * @brief GetMask Computes and returns the mask + * @return mitk::Image::Pointer of the generated mask + */ + mitk::Image::Pointer GetMask(); - void SetPlanarFigure(mitk::PlanarFigure::Pointer planarFigure); + void SetPlanarFigure(mitk::PlanarFigure::Pointer planarFigure); - /*We need an image because we need its geometry in order to give the points of the planar figure the correct world coordinates*/ - void SetImage(const mitk::Image::Pointer inputImage); + /** + * @brief SetImage: We need an image because we need its geometry in order to give the points of the planar figure the correct world coordinates + * @param inputImage + */ + void SetImage(const mitk::Image::Pointer inputImage); protected: - PlanarFigureMaskGenerator(){} + PlanarFigureMaskGenerator():Superclass(){} private: - void CalculateMask(); - - template < typename TPixel, unsigned int VImageDimension > - void InternalCalculateMaskFromPlanarFigure( - const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ); - - mitk::Image::Pointer extract2DImageSlice(unsigned int axis, unsigned int slice); - - bool GetPrincipalAxis(const BaseGeometry *geometry, Vector3D vector, - unsigned int &axis ); - - /** Connection from ITK to VTK */ - template - void ConnectPipelines(ITK_Exporter exporter, vtkSmartPointer importer) - { - importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); - - importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); - importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); - importer->SetSpacingCallback(exporter->GetSpacingCallback()); - importer->SetOriginCallback(exporter->GetOriginCallback()); - importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); - - importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); - - importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); - importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); - importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); - importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); - importer->SetCallbackUserData(exporter->GetCallbackUserData()); - } - - /** Connection from VTK to ITK */ - template - void ConnectPipelines(vtkSmartPointer exporter, ITK_Importer importer) - { - importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); - - importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); - importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); - importer->SetSpacingCallback(exporter->GetSpacingCallback()); - importer->SetOriginCallback(exporter->GetOriginCallback()); - importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); - - importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); - - importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); - importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); - importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); - importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); - importer->SetCallbackUserData(exporter->GetCallbackUserData()); - } - - - mitk::PlanarFigure::Pointer m_PlanarFigure; - typename itk::Image::Pointer m_InternalITKImageMask2D; - mitk::Image::Pointer m_InternalInputImage; - mitk::Image::Pointer m_InternalTimeSliceImage; - - unsigned int m_PlanarFigureAxis; + void CalculateMask(); + + template < typename TPixel, unsigned int VImageDimension > + void InternalCalculateMaskFromPlanarFigure( + const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ); + + mitk::Image::Pointer extract2DImageSlice(unsigned int axis, unsigned int slice); + + bool GetPrincipalAxis(const BaseGeometry *geometry, Vector3D vector, + unsigned int &axis ); + + /** Connection from ITK to VTK */ + template + void ConnectPipelines(ITK_Exporter exporter, vtkSmartPointer importer) + { + importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); + + importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); + importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); + importer->SetSpacingCallback(exporter->GetSpacingCallback()); + importer->SetOriginCallback(exporter->GetOriginCallback()); + importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); + + importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); + + importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); + importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); + importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); + importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); + importer->SetCallbackUserData(exporter->GetCallbackUserData()); + } + + /** Connection from VTK to ITK */ + template + void ConnectPipelines(vtkSmartPointer exporter, ITK_Importer importer) + { + importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); + + importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); + importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); + importer->SetSpacingCallback(exporter->GetSpacingCallback()); + importer->SetOriginCallback(exporter->GetOriginCallback()); + importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); + + importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); + + importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); + importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); + importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); + importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); + importer->SetCallbackUserData(exporter->GetCallbackUserData()); + } + + + mitk::PlanarFigure::Pointer m_PlanarFigure; + typename itk::Image::Pointer m_InternalITKImageMask2D; + mitk::Image::Pointer m_InternalInputImage; + mitk::Image::Pointer m_InternalTimeSliceImage; + + unsigned int m_PlanarFigureAxis; }; } #endif // MITKPLANARFIGUREMASKGENERATOR