diff --git a/Modules/Classification/CLMiniApps/CLPlanarFigureToNrrd.cpp b/Modules/Classification/CLMiniApps/CLPlanarFigureToNrrd.cpp index 488d45e581..da36d996ee 100644 --- a/Modules/Classification/CLMiniApps/CLPlanarFigureToNrrd.cpp +++ b/Modules/Classification/CLMiniApps/CLPlanarFigureToNrrd.cpp @@ -1,140 +1,140 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkCLPolyToNrrd_cpp #define mitkCLPolyToNrrd_cpp #include "time.h" #include #include #include "mitkCommandLineParser.h" #include #include #include #include #include typedef itk::Image< double, 3 > FloatImageType; typedef itk::Image< unsigned char, 3 > MaskImageType; struct MaskParameter { - mitk::Image::Pointer mask; + mitk::Image::ConstPointer mask; unsigned int axis; unsigned int slice; }; template < typename TPixel, unsigned int VImageDimension > void CreateNewMask(const itk::Image< TPixel, VImageDimension > *image, MaskParameter param, mitk::Image::Pointer &output) { int transform[3][2]; transform[0][0] = 1; transform[0][1] = 2; transform[1][0] = 0; transform[1][1] = 2; transform[2][0] = 0; transform[2][1] = 1; typedef itk::Image MaskType; typedef itk::Image Mask2DType; typename Mask2DType::Pointer mask = Mask2DType::New(); mitk::CastToItkImage(param.mask, mask); typename MaskType::Pointer mask3D = MaskType::New(); mask3D->SetRegions(image->GetLargestPossibleRegion()); mask3D->SetSpacing(image->GetSpacing()); mask3D->SetOrigin(image->GetOrigin()); mask3D->Allocate(); itk::ImageRegionIteratorWithIndex iter(mask3D, mask3D->GetLargestPossibleRegion()); while (!iter.IsAtEnd()) { auto index = iter.GetIndex(); iter.Set(0); if (index[param.axis] == param.slice) { Mask2DType::IndexType index2D; index2D[0] = index[transform[param.axis][0]]; index2D[1] = index[transform[param.axis][1]]; auto pixel = mask->GetPixel(index2D); iter.Set(pixel); } ++iter; } mitk::CastToMitkImage(mask3D, output); } int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); // required params parser.addArgument("planar", "p", mitkCommandLineParser::Directory, "Input Polydata", "Path to the input VTK polydata", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("image", "i", mitkCommandLineParser::Directory, "Input Image", "Image which defines the dimensions of the Segmentation", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.addArgument("output", "o", mitkCommandLineParser::File, "Output file", "Output file. ", us::Any(), false, false, false, mitkCommandLineParser::Input); // Miniapp Infos parser.setCategory("Classification Tools"); parser.setTitle("Planar Data to Nrrd Segmentation"); parser.setDescription("Creates a Nrrd segmentation based on a 2d-vtk polydata."); parser.setContributor("German Cancer Research Center (DKFZ)"); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) { return EXIT_FAILURE; } if ( parsedArgs.count("help") || parsedArgs.count("h")) { return EXIT_SUCCESS; } try { mitk::BaseData::Pointer data = mitk::IOUtil::Load(parsedArgs["planar"].ToString())[0]; mitk::PlanarFigure::Pointer planar = dynamic_cast(data.GetPointer()); mitk::Image::Pointer image = mitk::IOUtil::Load(parsedArgs["image"].ToString()); mitk::PlanarFigureMaskGenerator::Pointer pfMaskGen = mitk::PlanarFigureMaskGenerator::New(); pfMaskGen->SetPlanarFigure(planar); pfMaskGen->SetTimeStep(0); pfMaskGen->SetInputImage(image.GetPointer()); - mitk::Image::Pointer mask = pfMaskGen->GetMask(); + mitk::Image::ConstPointer mask = pfMaskGen->GetMask(); unsigned int axis = pfMaskGen->GetPlanarFigureAxis(); unsigned int slice = pfMaskGen->GetPlanarFigureSlice(); //itk::Image::IndexType index; mitk::Image::Pointer fullMask; MaskParameter param; param.slice = slice; param.axis = axis; param.mask = mask; AccessByItk_2(image, CreateNewMask, param, fullMask); std::string saveAs = parsedArgs["output"].ToString(); MITK_INFO << "Save as: " << saveAs; - mitk::IOUtil::Save(pfMaskGen->GetMask(), saveAs); + mitk::IOUtil::Save(mask, saveAs); mitk::IOUtil::Save(fullMask, saveAs); return 0; } catch (...) { return EXIT_FAILURE; } } #endif diff --git a/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp b/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp index b993b0f89e..fff3c3eee2 100644 --- a/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp +++ b/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp @@ -1,649 +1,650 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "mitkImageStatisticsCalculator.h" #include "itkMultiGaussianImageSource.h" #include "mitkTestingMacros.h" #include "mitkImageCast.h" #include #include #include #include #include #include #include /** \section hotspotCalculationTestCases Testcases To see the different Hotspot-Testcases have a look at the \ref hotspottestdoc. Note from an intensive session of checking the test results: - itk::MultiGaussianImageSource needs a review - the test idea is ok, but the combination of XML files for parameters and MultiGaussianImageSource has serious flaws - the XML file should contain exactly the parameters that MultiGaussianImageSource requires - in contrast, now the XML file mentions index coordinates for gaussian centers while the MultiGaussianImageSource expects world coordinates - this requires a transformation (index * spacing assuming no rotation) that was actually broken until recently */ struct mitkImageStatisticsHotspotTestClass { /** \brief Test parameters for one test case. Describes all aspects of a single test case: - parameters to generate a test image - parameters of a ROI that describes where to calculate statistics - expected statistics results */ struct Parameters { public: // XML-Tag /** \brief XML-Tag "image-rows": size of x-dimension */ int m_ImageRows; /** \brief XML-Tag "image-columns": size of y-dimension */ int m_ImageColumns; /** \brief XML-Tag "image-slices": size of z-dimension */ int m_ImageSlices; /** \brief XML-Tag "numberOfGaussians": number of used gauss-functions */ int m_NumberOfGaussian; /** \brief XML-Tags "spacingX", "spacingY", "spacingZ": spacing of image in every direction */ double m_Spacing[3]; /** \brief XML-Tag "entireHotSpotInImage" */ unsigned int m_EntireHotspotInImage; // XML-Tag /** \brief XML-Tag "centerIndexX: gaussian parameter \warning This parameter READS the centerIndexX parameter from file and is THEN MISUSED to calculate some position in world coordinates, so we require double. */ std::vector m_CenterX; /** \brief XML-Tag "centerIndexY: gaussian parameter \warning This parameter READS the centerIndexX parameter from file and is THEN MISUSED to calculate some position in world coordinates, so we require double. */ std::vector m_CenterY; /** \brief XML-Tag "centerIndexZ: gaussian parameter \warning This parameter READS the centerIndexX parameter from file and is THEN MISUSED to calculate some position in world coordinates, so we require double. */ std::vector m_CenterZ; /** \brief XML-Tag "deviationX: gaussian parameter */ std::vector m_SigmaX; /** \brief XML-Tag "deviationY: gaussian parameter */ std::vector m_SigmaY; /** \brief XML-Tag "deviationZ: gaussian parameter */ std::vector m_SigmaZ; /** \brief XML-Tag "altitude: gaussian parameter */ std::vector m_Altitude; // XML-Tag /** \brief XML-Tag "numberOfLabels": number of different labels which appear in the mask */ unsigned int m_NumberOfLabels; /** \brief XML-Tag "hotspotRadiusInMM": radius of hotspot */ double m_HotspotRadiusInMM; // XML-Tag /** \brief XML-Tag "maximumSizeX": maximum position of ROI in x-dimension */ vnl_vector m_MaxIndexX; /** \brief XML-Tag "minimumSizeX": minimum position of ROI in x-dimension */ vnl_vector m_MinIndexX; /** \brief XML-Tag "maximumSizeX": maximum position of ROI in y-dimension */ vnl_vector m_MaxIndexY; /** \brief XML-Tag "minimumSizeX": minimum position of ROI in y-dimension */ vnl_vector m_MinIndexY; /** \brief XML-Tag "maximumSizeX": maximum position of ROI in z-dimension */ vnl_vector m_MaxIndexZ; /** \brief XML-Tag "minimumSizeX": minimum position of ROI in z-dimension */ vnl_vector m_MinIndexZ; /** \brief XML-Tag "label": value of label */ vnl_vector m_Label; //XML-Tag /** \brief XML-Tag "minimum": minimum inside hotspot */ vnl_vector m_HotspotMin; /** \brief XML-Tag "maximum": maximum inside hotspot */ vnl_vector m_HotspotMax; /** \brief XML-Tag "mean": mean value of hotspot */ vnl_vector m_HotspotMean; /** \brief XML-Tag "maximumIndexX": x-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMaxIndexX; /** \brief XML-Tag "maximumIndexX": y-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMaxIndexY; /** \brief XML-Tag "maximumIndexX": z-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMaxIndexZ; /** \brief XML-Tag "maximumIndexX": x-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMinIndexX; /** \brief XML-Tag "maximumIndexX": y-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMinIndexY; /** \brief XML-Tag "maximumIndexX": z-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMinIndexZ; /** \brief XML-Tag "maximumIndexX": x-coordinate of hotspot-location */ vnl_vector m_HotspotIndexX; /** \brief XML-Tag "maximumIndexX": y-coordinate of hotspot-location */ vnl_vector m_HotspotIndexY; /** \brief XML-Tag "maximumIndexX": z-coordinate of hotspot-location */ vnl_vector m_HotspotIndexZ; }; /** \brief Find/Convert integer attribute in itk::DOMNode. */ static int GetIntegerAttribute(itk::DOMNode* domNode, const std::string& tag) { assert(domNode); MITK_TEST_CONDITION_REQUIRED( domNode->HasAttribute(tag), "Tag '" << tag << "' is defined in test parameters" ); std::string attributeValue = domNode->GetAttribute(tag); int resultValue; try { //MITK_TEST_OUTPUT( << "Converting tag value '" << attributeValue << "' for tag '" << tag << "' to integer"); std::stringstream(attributeValue) >> resultValue; return resultValue; } catch(std::exception& /*e*/) { MITK_TEST_CONDITION_REQUIRED(false, "Convert tag value '" << attributeValue << "' for tag '" << tag << "' to integer"); return 0; // just to satisfy compiler } } /** \brief Find/Convert double attribute in itk::DOMNode. */ static double GetDoubleAttribute(itk::DOMNode* domNode, const std::string& tag) { assert(domNode); MITK_TEST_CONDITION_REQUIRED( domNode->HasAttribute(tag), "Tag '" << tag << "' is defined in test parameters" ); std::string attributeValue = domNode->GetAttribute(tag); double resultValue; try { //MITK_TEST_OUTPUT( << "Converting tag value '" << attributeValue << "' for tag '" << tag << "' to double"); std::stringstream(attributeValue) >> resultValue; return resultValue; } catch(std::exception& /*e*/) { MITK_TEST_CONDITION_REQUIRED(false, "Convert tag value '" << attributeValue << "' for tag '" << tag << "' to double"); return 0.0; // just to satisfy compiler } } /** \brief Read XML file describing the test parameters. Reads XML file given in first commandline parameter in order to construct a Parameters structure. The XML file should be structurs as the following example, i.e. we describe the three test aspects of Parameters in four different tags, with all the details described as tag attributes. */ /** \verbatim \endverbatim */ static Parameters ParseParameters(int argc, char* argv[]) { MITK_TEST_CONDITION_REQUIRED(argc == 2, "Test is invoked with exactly 1 parameter (XML parameters file)"); MITK_INFO << "Reading parameters from file '" << argv[1] << "'"; std::string filename = argv[1]; Parameters result; itk::DOMNodeXMLReader::Pointer xmlReader = itk::DOMNodeXMLReader::New(); xmlReader->SetFileName( filename ); try { xmlReader->Update(); itk::DOMNode::Pointer domRoot = xmlReader->GetOutput(); typedef std::vector NodeList; NodeList testimages; domRoot->GetChildren("testimage", testimages); MITK_TEST_CONDITION_REQUIRED( testimages.size() == 1, "One test image defined" ) itk::DOMNode* testimage = testimages[0]; result.m_ImageRows = GetIntegerAttribute( testimage, "image-rows" ); result.m_ImageColumns = GetIntegerAttribute( testimage, "image-columns" ); result.m_ImageSlices = GetIntegerAttribute( testimage, "image-slices" ); result.m_NumberOfGaussian = GetIntegerAttribute( testimage, "numberOfGaussians" ); result.m_Spacing[0] = GetDoubleAttribute(testimage, "spacingX"); result.m_Spacing[1] = GetDoubleAttribute(testimage, "spacingY"); result.m_Spacing[2] = GetDoubleAttribute(testimage, "spacingZ"); result.m_EntireHotspotInImage = GetIntegerAttribute( testimage, "entireHotSpotInImage" ); MITK_TEST_OUTPUT( << "Read size parameters (x,y,z): " << result.m_ImageRows << "," << result.m_ImageColumns << "," << result.m_ImageSlices); MITK_TEST_OUTPUT( << "Read spacing parameters (x,y,z): " << result.m_Spacing[0] << "," << result.m_Spacing[1] << "," << result.m_Spacing[2]); NodeList gaussians; testimage->GetChildren("gaussian", gaussians); MITK_TEST_CONDITION_REQUIRED( gaussians.size() >= 1, "At least one gaussian is defined" ) result.m_CenterX.resize(result.m_NumberOfGaussian); result.m_CenterY.resize(result.m_NumberOfGaussian); result.m_CenterZ.resize(result.m_NumberOfGaussian); result.m_SigmaX.resize(result.m_NumberOfGaussian); result.m_SigmaY.resize(result.m_NumberOfGaussian); result.m_SigmaZ.resize(result.m_NumberOfGaussian); result.m_Altitude.resize(result.m_NumberOfGaussian); for(int i = 0; i < result.m_NumberOfGaussian ; ++i) { itk::DOMNode* gaussian = gaussians[i]; result.m_CenterX[i] = GetIntegerAttribute(gaussian, "centerIndexX"); result.m_CenterY[i] = GetIntegerAttribute(gaussian, "centerIndexY"); result.m_CenterZ[i] = GetIntegerAttribute(gaussian, "centerIndexZ"); result.m_SigmaX[i] = GetDoubleAttribute(gaussian, "deviationX"); result.m_SigmaY[i] = GetDoubleAttribute(gaussian, "deviationY"); result.m_SigmaZ[i] = GetDoubleAttribute(gaussian, "deviationZ"); result.m_Altitude[i] = GetDoubleAttribute(gaussian, "altitude"); result.m_CenterX[i] = result.m_CenterX[i] * result.m_Spacing[0]; result.m_CenterY[i] = result.m_CenterY[i] * result.m_Spacing[1]; result.m_CenterZ[i] = result.m_CenterZ[i] * result.m_Spacing[2]; result.m_SigmaX[i] = result.m_SigmaX[i] * result.m_Spacing[0]; result.m_SigmaY[i] = result.m_SigmaY[i] * result.m_Spacing[1]; result.m_SigmaZ[i] = result.m_SigmaZ[i] * result.m_Spacing[2]; } NodeList segmentations; domRoot->GetChildren("segmentation", segmentations); MITK_TEST_CONDITION_REQUIRED( segmentations.size() == 1, "One segmentation defined"); itk::DOMNode* segmentation = segmentations[0]; result.m_NumberOfLabels = GetIntegerAttribute(segmentation, "numberOfLabels"); result.m_HotspotRadiusInMM = GetDoubleAttribute(segmentation, "hotspotRadiusInMM"); // read ROI parameters, fill result structure NodeList rois; segmentation->GetChildren("roi", rois); MITK_TEST_CONDITION_REQUIRED( rois.size() >= 1, "At least one ROI defined" ) result.m_MaxIndexX.set_size(result.m_NumberOfLabels); result.m_MinIndexX.set_size(result.m_NumberOfLabels); result.m_MaxIndexY.set_size(result.m_NumberOfLabels); result.m_MinIndexY.set_size(result.m_NumberOfLabels); result.m_MaxIndexZ.set_size(result.m_NumberOfLabels); result.m_MinIndexZ.set_size(result.m_NumberOfLabels); result.m_Label.set_size(result.m_NumberOfLabels); for(unsigned int i = 0; i < rois.size(); ++i) { result.m_MaxIndexX[i] = GetIntegerAttribute(rois[i], "maximumIndexX"); result.m_MinIndexX[i] = GetIntegerAttribute(rois[i], "minimumIndexX"); result.m_MaxIndexY[i] = GetIntegerAttribute(rois[i], "maximumIndexY"); result.m_MinIndexY[i] = GetIntegerAttribute(rois[i], "minimumIndexY"); result.m_MaxIndexZ[i] = GetIntegerAttribute(rois[i], "maximumIndexZ"); result.m_MinIndexZ[i] = GetIntegerAttribute(rois[i], "minimumIndexZ"); result.m_Label[i] = GetIntegerAttribute(rois[i], "label"); } // read statistic parameters, fill result structure NodeList statistics; domRoot->GetChildren("statistic", statistics); MITK_TEST_CONDITION_REQUIRED( statistics.size() >= 1 , "At least one statistic defined" ) MITK_TEST_CONDITION_REQUIRED( statistics.size() == rois.size(), "Same number of rois and corresponding statistics defined"); result.m_HotspotMin.set_size(statistics.size()); result.m_HotspotMax.set_size(statistics.size()); result.m_HotspotMean.set_size(statistics.size()); result.m_HotspotMinIndexX.set_size(statistics.size()); result.m_HotspotMinIndexY.set_size(statistics.size()); result.m_HotspotMinIndexZ.set_size(statistics.size()); result.m_HotspotMaxIndexX.set_size(statistics.size()); result.m_HotspotMaxIndexY.set_size(statistics.size()); result.m_HotspotMaxIndexZ.set_size(statistics.size()); result.m_HotspotIndexX.set_size(statistics.size()); result.m_HotspotIndexY.set_size(statistics.size()); result.m_HotspotIndexZ.set_size(statistics.size()); for(unsigned int i = 0; i < statistics.size(); ++i) { result.m_HotspotMin[i] = GetDoubleAttribute(statistics[i], "minimum"); result.m_HotspotMax[i] = GetDoubleAttribute(statistics[i], "maximum"); result.m_HotspotMean[i] = GetDoubleAttribute(statistics[i], "mean"); result.m_HotspotMinIndexX[i] = GetIntegerAttribute(statistics[i], "minimumIndexX"); result.m_HotspotMinIndexY[i] = GetIntegerAttribute(statistics[i], "minimumIndexY"); result.m_HotspotMinIndexZ[i] = GetIntegerAttribute(statistics[i], "minimumIndexZ"); result.m_HotspotMaxIndexX[i] = GetIntegerAttribute(statistics[i], "maximumIndexX"); result.m_HotspotMaxIndexY[i] = GetIntegerAttribute(statistics[i], "maximumIndexY"); result.m_HotspotMaxIndexZ[i] = GetIntegerAttribute(statistics[i], "maximumIndexZ"); result.m_HotspotIndexX[i] = GetIntegerAttribute(statistics[i], "hotspotIndexX"); result.m_HotspotIndexY[i] = GetIntegerAttribute(statistics[i], "hotspotIndexY"); result.m_HotspotIndexZ[i] = GetIntegerAttribute(statistics[i], "hotspotIndexZ"); } } catch (std::exception& e) { MITK_TEST_CONDITION_REQUIRED(false, "Reading test parameters from XML file. Error message: " << e.what()); } return result; } /** \brief Generate an image that contains a couple of 3D gaussian distributions. Uses the given parameters to produce a test image using class MultiGaussianImageSource. */ static mitk::Image::Pointer BuildTestImage(const Parameters& testParameters) { mitk::Image::Pointer result; typedef double PixelType; const int Dimension = 3; typedef itk::Image ImageType; ImageType::Pointer image = ImageType::New(); typedef itk::MultiGaussianImageSource< ImageType > MultiGaussianImageSource; MultiGaussianImageSource::Pointer gaussianGenerator = MultiGaussianImageSource::New(); ImageType::SizeValueType size[3]; size[0] = testParameters.m_ImageColumns; size[1] = testParameters.m_ImageRows; size[2] = testParameters.m_ImageSlices; itk::MultiGaussianImageSource::VectorType centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec; for(int i = 0; i < testParameters.m_NumberOfGaussian; ++i) { centerXVec.push_back(testParameters.m_CenterX[i]); centerYVec.push_back(testParameters.m_CenterY[i]); centerZVec.push_back(testParameters.m_CenterZ[i]); sigmaXVec.push_back(testParameters.m_SigmaX[i]); sigmaYVec.push_back(testParameters.m_SigmaY[i]); sigmaZVec.push_back(testParameters.m_SigmaZ[i]); altitudeVec.push_back(testParameters.m_Altitude[i]); } ImageType::SpacingType spacing; for( int i = 0; i < Dimension; ++i ) spacing[i] = testParameters.m_Spacing[i]; gaussianGenerator->SetSize( size ); gaussianGenerator->SetSpacing( spacing ); gaussianGenerator->SetRadius(testParameters.m_HotspotRadiusInMM); gaussianGenerator->SetNumberOfGausssians(testParameters.m_NumberOfGaussian); gaussianGenerator->AddGaussian(centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec); gaussianGenerator->Update(); image = gaussianGenerator->GetOutput(); mitk::CastToMitkImage(image, result); return result; } /** \brief Calculates hotspot statistics for given test image and ROI parameters. Uses ImageStatisticsCalculator to find a hotspot in a defined ROI within the given image. */ static mitk::ImageStatisticsContainer::ImageStatisticsObject CalculateStatistics(mitk::Image* image, const Parameters& testParameters, unsigned int label) { const unsigned int Dimension = 3; typedef itk::Image MaskImageType; MaskImageType::Pointer mask = MaskImageType::New(); MaskImageType::SizeType size; MaskImageType::SpacingType spacing; MaskImageType::IndexType start; mitk::ImageStatisticsCalculator::Pointer statisticsCalculator = mitk::ImageStatisticsCalculator::New(); statisticsCalculator->SetInputImage(image); mitk::Image::Pointer mitkMaskImage; if((testParameters.m_MaxIndexX[label] > testParameters.m_MinIndexX[label] && testParameters.m_MinIndexX[label] >= 0) && (testParameters.m_MaxIndexY[label] > testParameters.m_MinIndexY[label] && testParameters.m_MinIndexY[label] >= 0) && (testParameters.m_MaxIndexZ[label] > testParameters.m_MinIndexZ[label] && testParameters.m_MinIndexZ[label] >= 0)) { for(unsigned int i = 0; i < Dimension; ++i) { start[i] = 0; spacing[i] = testParameters.m_Spacing[i]; } size[0] = testParameters.m_ImageColumns; size[1] = testParameters.m_ImageRows; size[2] = testParameters.m_ImageSlices; MaskImageType::RegionType region; region.SetIndex(start); region.SetSize(size); mask->SetSpacing(spacing); mask->SetRegions(region); mask->Allocate(); typedef itk::ImageRegionIteratorWithIndex MaskImageIteratorType; MaskImageIteratorType maskIt(mask, region); for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { maskIt.Set(0); } for(unsigned int i = 0; i < testParameters.m_NumberOfLabels; ++i) { for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { MaskImageType::IndexType index = maskIt.GetIndex(); if((index[0] >= testParameters.m_MinIndexX[i] && index[0] <= testParameters.m_MaxIndexX[i] ) && (index[1] >= testParameters.m_MinIndexY[i] && index[1] <= testParameters.m_MaxIndexY[i] ) && (index[2] >= testParameters.m_MinIndexZ[i] && index[2] <= testParameters.m_MaxIndexZ[i] )) { maskIt.Set(testParameters.m_Label[i]); } } } mitk::CastToMitkImage(mask, mitkMaskImage); mitk::ImageMaskGenerator::Pointer imgMaskGen = mitk::ImageMaskGenerator::New(); + imgMaskGen->SetInputImage(image); imgMaskGen->SetImageMask(mitkMaskImage); mitk::HotspotMaskGenerator::Pointer hotspotMaskGen = mitk::HotspotMaskGenerator::New(); hotspotMaskGen->SetInputImage(image); hotspotMaskGen->SetLabel(testParameters.m_Label[label]); hotspotMaskGen->SetMask(imgMaskGen.GetPointer()); hotspotMaskGen->SetHotspotRadiusInMM(testParameters.m_HotspotRadiusInMM); if(testParameters.m_EntireHotspotInImage == 1) { MITK_INFO << "Hotspot must be completly inside image"; hotspotMaskGen->SetHotspotMustBeCompletelyInsideImage(true); } else { MITK_INFO << "Hotspot must not be completly inside image"; hotspotMaskGen->SetHotspotMustBeCompletelyInsideImage(false); } statisticsCalculator->SetMask(hotspotMaskGen.GetPointer()); MITK_DEBUG << "Masking is set to hotspot+image mask"; } else { mitk::HotspotMaskGenerator::Pointer hotspotMaskGen = mitk::HotspotMaskGenerator::New(); hotspotMaskGen->SetInputImage(image); hotspotMaskGen->SetHotspotRadiusInMM(testParameters.m_HotspotRadiusInMM); if(testParameters.m_EntireHotspotInImage == 1) { MITK_INFO << "Hotspot must be completly inside image"; hotspotMaskGen->SetHotspotMustBeCompletelyInsideImage(true); } else { MITK_INFO << "Hotspot must not be completly inside image"; hotspotMaskGen->SetHotspotMustBeCompletelyInsideImage(false); } MITK_DEBUG << "Masking is set to hotspot only"; } return statisticsCalculator->GetStatistics()->GetStatisticsForTimeStep(0); } static void ValidateStatisticsItem(const std::string& label, double testvalue, double reference, double tolerance) { double diff = ::fabs(reference - testvalue); MITK_TEST_CONDITION( diff < tolerance, "'" << label << "' value close enough to reference value " "(value=" << testvalue << ", reference=" << reference << ", diff=" << diff << ")" ); } static void ValidateStatisticsItem(const std::string& label, const vnl_vector& testvalue, const vnl_vector& reference) { double diffX = ::fabs(double(testvalue[0] - reference[0])); double diffY = ::fabs(double(testvalue[1] - reference[1])); double diffZ = ::fabs(double(testvalue[2] - reference[2])); std::stringstream testPosition; testPosition << testvalue[0] << "," << testvalue[1] << "," << testvalue[2]; std::stringstream referencePosition; referencePosition << reference[0] << "," << reference[1] << "," << reference[2]; MITK_TEST_CONDITION( diffX < mitk::eps && diffY < mitk::eps && diffZ < mitk::eps, "'" << label << "' close enough to reference value " << "(value=[" << testPosition.str() << "]," << " reference=[" << referencePosition.str() << "]"); } /** \brief Compares calculated against actual statistics values. Checks validness of all statistics aspects. Lets test fail if any aspect is not sufficiently equal. */ static void ValidateStatistics(const mitk::ImageStatisticsContainer::ImageStatisticsObject hotspotStatistics, const Parameters& testParameters, unsigned int label) { // check all expected test result against actual results double eps = 0.25; // value above the largest tested difference auto mean = hotspotStatistics.GetValueConverted(mitk::ImageStatisticsConstants::MEAN()); auto max = hotspotStatistics.GetValueConverted(mitk::ImageStatisticsConstants::MAXIMUM()); auto min = hotspotStatistics.GetValueConverted(mitk::ImageStatisticsConstants::MINIMUM()); ValidateStatisticsItem("Hotspot mean", mean, testParameters.m_HotspotMean[label], eps); ValidateStatisticsItem("Hotspot maximum", max, testParameters.m_HotspotMax[label], eps); ValidateStatisticsItem("Hotspot minimum", min, testParameters.m_HotspotMin[label], eps); vnl_vector referenceHotspotCenterIndex; referenceHotspotCenterIndex.set_size(3); referenceHotspotCenterIndex[0] = testParameters.m_HotspotIndexX[label]; referenceHotspotCenterIndex[1] = testParameters.m_HotspotIndexY[label]; referenceHotspotCenterIndex[2] = testParameters.m_HotspotIndexZ[label]; // ValidateStatisticsItem("Hotspot center position", statistics.GetHotspotStatistics().GetHotspotIndex(), referenceHotspotCenterIndex); TODO: new image statistics calculator does not give hotspot position // TODO we do not test minimum/maximum positions within the peak/hotspot region, because // these positions are not unique, i.e. there are multiple valid minima/maxima positions. // One solution would be to modify the test cases in order to achive clear positions. // The BETTER/CORRECT solution would be to change the singular position into a set of positions / a region } }; /** \brief Verifies that hotspot statistics part of ImageStatisticsCalculator. The test reads parameters from an XML-file to generate a test-image, calculates the hotspot statistics of the image and checks if the calculated statistics are the same as the specified values of the XML-file. */ int mitkImageStatisticsHotspotTest(int argc, char* argv[]) { MITK_TEST_BEGIN("mitkImageStatisticsHotspotTest") try { mitkImageStatisticsHotspotTestClass::Parameters parameters = mitkImageStatisticsHotspotTestClass::ParseParameters(argc,argv); mitk::Image::Pointer image = mitkImageStatisticsHotspotTestClass::BuildTestImage(parameters); MITK_TEST_CONDITION_REQUIRED( image.IsNotNull(), "Generate test image" ); for(unsigned int label = 0; label < parameters.m_NumberOfLabels; ++label) { mitk::ImageStatisticsContainer::ImageStatisticsObject statistics = mitkImageStatisticsHotspotTestClass::CalculateStatistics(image, parameters, label); mitkImageStatisticsHotspotTestClass::ValidateStatistics(statistics, parameters, label); std::cout << std::endl; } } catch (std::exception& e) { std::cout << "Error: " << e.what() << std::endl; MITK_TEST_CONDITION_REQUIRED( false, "Exception occurred during test execution: " << e.what() ); } catch(...) { MITK_TEST_CONDITION_REQUIRED( false, "Exception occurred during test execution." ); } MITK_TEST_END() } diff --git a/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp b/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp index 30c805a19c..35d4a30a8e 100644 --- a/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp @@ -1,612 +1,618 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include #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_TimeStep = 0; m_InternalMask = mitk::Image::New(); m_InternalMaskUpdateTime = 0; } void HotspotMaskGenerator::SetInputImage(mitk::Image::Pointer inputImage) { if (inputImage != m_inputImage) { m_inputImage = inputImage; m_ConvolutionImageMaxIndex.set_size(inputImage->GetDimension()); m_ConvolutionImageMinIndex.set_size(inputImage->GetDimension()); this->Modified(); } } void HotspotMaskGenerator::SetMask(MaskGenerator::Pointer mask) { if (mask != m_Mask) { m_Mask = mask; this->Modified(); } } HotspotMaskGenerator::~HotspotMaskGenerator() { } void HotspotMaskGenerator::SetHotspotRadiusInMM(double radiusInMillimeter) { if(radiusInMillimeter != m_HotspotRadiusinMM) { m_HotspotRadiusinMM = radiusInMillimeter; this->Modified(); } } const double& HotspotMaskGenerator::GetHotspotRadiusinMM() const { return m_HotspotRadiusinMM; } bool HotspotMaskGenerator::GetHotspotMustBeCompletelyInsideImage() const { return m_HotspotMustBeCompletelyInsideImage; } void HotspotMaskGenerator::SetHotspotMustBeCompletelyInsideImage(bool mustBeCompletelyInImage) { if (m_HotspotMustBeCompletelyInsideImage != mustBeCompletelyInImage) { m_HotspotMustBeCompletelyInsideImage = mustBeCompletelyInImage; this->Modified(); } } - mitk::Image::Pointer HotspotMaskGenerator::GetMask() + mitk::Image::ConstPointer HotspotMaskGenerator::GetMask() { if (IsUpdateRequired()) { if ( m_inputImage.IsNull() ) { throw std::runtime_error( "Error: image empty!" ); } if ( m_TimeStep >= m_inputImage->GetTimeSteps() ) { throw std::runtime_error( "Error: invalid time step!" ); } mitk::ImageTimeSelector::Pointer imageTimeSelector = mitk::ImageTimeSelector::New(); imageTimeSelector->SetInput( m_inputImage ); 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(); + mitk::Image::ConstPointer timeSliceMask = m_Mask->GetMask(); if ( m_internalImage->GetDimension() == 3 ) { - CastToItkImage(timeSliceMask, m_internalMask3D); - AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D, m_Label); + itk::Image::Pointer noneConstMaskImage; //needed to work arround the fact that CastToItkImage currently does not support const itk images. + CastToItkImage(timeSliceMask, noneConstMaskImage); + m_internalMask3D = noneConstMaskImage; + AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D.GetPointer(), m_Label); } else if ( m_internalImage->GetDimension() == 2 ) { - CastToItkImage(timeSliceMask, m_internalMask2D); - AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D, m_Label); + itk::Image::Pointer noneConstMaskImage; //needed to work arround the fact that CastToItkImage currently does not support const itk images. + CastToItkImage(timeSliceMask, noneConstMaskImage); + m_internalMask2D = noneConstMaskImage; + AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D.GetPointer(), 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); + AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D.GetPointer(), m_Label); } else if ( m_internalImage->GetDimension() == 2 ) { - AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D, m_Label); + AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D.GetPointer(), m_Label); } else { throw std::runtime_error( "Error: invalid image dimension" ); } } this->Modified(); } m_InternalMaskUpdateTime = m_InternalMask->GetMTime(); return m_InternalMask; } void HotspotMaskGenerator::SetTimeStep(unsigned int timeStep) { if (m_TimeStep != timeStep) { m_TimeStep = timeStep; } } void HotspotMaskGenerator::SetLabel(unsigned short label) { if (label != m_Label) { m_Label = label; this->Modified(); } } vnl_vector HotspotMaskGenerator::GetConvolutionImageMinIndex() { this->GetMask(); // make sure we are up to date return m_ConvolutionImageMinIndex; } vnl_vector HotspotMaskGenerator::GetHotspotIndex() { this->GetMask(); // make sure we are up to date return m_ConvolutionImageMaxIndex; } template HotspotMaskGenerator::ImageExtrema HotspotMaskGenerator::CalculateExtremaWorld( const itk::Image* inputImage, - typename itk::Image::Pointer maskImage, + const itk::Image* 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) { itk::IndexValueType 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::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; 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]; } double 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, + const itk::Image* maskImage, unsigned int label) { typedef itk::Image< TPixel, VImageDimension > InputImageType; typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; 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()"); } + typename MaskImageType::ConstPointer usedMask = maskImage; // 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(); + auto defaultMask = 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); + defaultMask->SetRegions(maskRegion); + defaultMask->Allocate(); + defaultMask->SetOrigin(maskOrigin); + defaultMask->SetSpacing(maskSpacing); + defaultMask->SetDirection(maskDirection); - maskImage->FillBuffer(1); + defaultMask->FillBuffer(1); + usedMask = defaultMask; 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); + ImageExtrema convolutionImageInformation = CalculateExtremaWorld(convolutionImage.GetPointer(), usedMask.GetPointer(), 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 = 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; m_ConvolutionImageMaxIndex = convolutionImageInformation.MaxIndex; m_ConvolutionImageMinIndex = convolutionImageInformation.MinIndex; } } bool HotspotMaskGenerator::IsUpdateRequired() const { unsigned long thisClassTimeStamp = this->GetMTime(); unsigned long internalMaskTimeStamp = m_InternalMask->GetMTime(); unsigned long maskGeneratorTimeStamp = m_Mask->GetMTime(); unsigned long inputImageTimeStamp = m_inputImage->GetMTime(); if (thisClassTimeStamp > m_InternalMaskUpdateTime) // inputs have changed { return true; } if (m_InternalMaskUpdateTime < maskGeneratorTimeStamp || m_InternalMaskUpdateTime < inputImageTimeStamp) // mask image has changed outside of this class { return true; } if (internalMaskTimeStamp > m_InternalMaskUpdateTime) // internal mask has been changed outside of this class { return true; } return false; } } diff --git a/Modules/ImageStatistics/mitkHotspotMaskGenerator.h b/Modules/ImageStatistics/mitkHotspotMaskGenerator.h index 724538410e..e23299db6a 100644 --- a/Modules/ImageStatistics/mitkHotspotMaskGenerator.h +++ b/Modules/ImageStatistics/mitkHotspotMaskGenerator.h @@ -1,179 +1,180 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef MITKHOTSPOTCALCULATOR_H #define MITKHOTSPOTCALCULATOR_H #include #include #include #include #include #include #include #include namespace mitk { /** * @brief The HotspotMaskGenerator class is used when a hotspot has to be found in an image. A hotspot is * the region of the image where the mean intensity is maximal (=brightest spot). It is usually used in PET scans. * The identification of the hotspot is done as follows: First a cubic (or circular, if image is 2d) * mask of predefined size is generated. This mask is then convolved with the input image (in fourier domain). * The maximum value of the convolved image then corresponds to the hotspot. * If a maskGenerator is set, only the pixels of the convolved image where the corresponding mask is == @a label * are searched for the maximum value. */ class MITKIMAGESTATISTICS_EXPORT HotspotMaskGenerator: public MaskGenerator { public: /** Standard Self typedef */ typedef HotspotMaskGenerator 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(HotspotMaskGenerator, MaskGenerator); /** @brief Set the input image. Required for this class */ void SetInputImage(mitk::Image::Pointer inputImage); /** @brief Set a mask (can be nullptr if no mask is desired) */ void SetMask(MaskGenerator::Pointer mask); /** @brief Set the radius of the hotspot (in MM) */ void SetHotspotRadiusInMM(double radiusInMillimeter); const double& GetHotspotRadiusinMM() const; /** @brief Define whether the hotspot must be completely inside the image. Default is true */ void SetHotspotMustBeCompletelyInsideImage(bool hotspotCompletelyInsideImage); bool GetHotspotMustBeCompletelyInsideImage() const; /** @brief If a maskGenerator is set, this detemines which mask value is used */ void SetLabel(unsigned short label); /** @brief Computes and returns the hotspot mask. The hotspot mask has the same size as the input image. The hopspot has value 1, the remaining pixels are set to 0 */ - mitk::Image::Pointer GetMask() override; + mitk::Image::ConstPointer GetMask() override; /** @brief Returns the image index where the hotspot is located */ vnl_vector GetHotspotIndex(); /** @brief Returns the index where the convolution image is minimal (darkest spot in image) */ vnl_vector GetConvolutionImageMinIndex(); /** * @brief SetTimeStep is used to set the time step for which the mask is to be generated * @param timeStep */ void SetTimeStep(unsigned int timeStep) override; protected: HotspotMaskGenerator(); ~HotspotMaskGenerator() override; class ImageExtrema { public: bool Defined; double Max; double Min; vnl_vector MaxIndex; vnl_vector MinIndex; ImageExtrema() :Defined(false) ,Max(itk::NumericTraits::min()) ,Min(itk::NumericTraits::max()) { } }; private: /** \brief Returns size of convolution kernel depending on spacing and radius. */ template itk::Size CalculateConvolutionKernelSize(double spacing[VImageDimension], double radiusInMM); /** \brief Generates image of kernel which is needed for convolution. */ template itk::SmartPointer< itk::Image > GenerateHotspotSearchConvolutionKernel(double spacing[VImageDimension], double radiusInMM); /** \brief Convolves image with spherical kernel image. Used for hotspot calculation. */ template itk::SmartPointer< itk::Image > GenerateConvolutionImage( const itk::Image* inputImage ); /** \brief Fills pixels of the spherical hotspot mask. */ template < typename TPixel, unsigned int VImageDimension> void FillHotspotMaskPixels( itk::Image* maskImage, itk::Point sphereCenter, double sphereRadiusInMM); /** \brief */ template void CalculateHotspotMask(itk::Image* inputImage, - typename itk::Image::Pointer maskImage, + const itk::Image* maskImage, unsigned int label); template ImageExtrema CalculateExtremaWorld( const itk::Image* inputImage, - typename itk::Image::Pointer maskImage, + const itk::Image* maskImage, double neccessaryDistanceToImageBorderInMM, unsigned int label); bool IsUpdateRequired() const; HotspotMaskGenerator(const HotspotMaskGenerator &); HotspotMaskGenerator & operator=(const HotspotMaskGenerator &); MaskGenerator::Pointer m_Mask; + mitk::Image::Pointer m_InternalMask; mitk::Image::Pointer m_internalImage; - itk::Image::Pointer m_internalMask2D; - itk::Image::Pointer m_internalMask3D; + itk::Image::ConstPointer m_internalMask2D; + itk::Image::ConstPointer m_internalMask3D; double m_HotspotRadiusinMM; bool m_HotspotMustBeCompletelyInsideImage; unsigned short m_Label; vnl_vector m_ConvolutionImageMinIndex, m_ConvolutionImageMaxIndex; unsigned long m_InternalMaskUpdateTime; }; } #endif // MITKHOTSPOTCALCULATOR diff --git a/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.cpp b/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.cpp index 3f35368492..ad1787bc84 100644 --- a/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.cpp @@ -1,132 +1,132 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include #include #include #include #include #include namespace mitk { void IgnorePixelMaskGenerator::SetIgnoredPixelValue(RealType pixelValue) { if (pixelValue != m_IgnoredPixelValue) { m_IgnoredPixelValue = pixelValue; this->Modified(); } } void IgnorePixelMaskGenerator::SetTimeStep(unsigned int timeStep) { if (m_TimeStep != timeStep) { m_TimeStep = timeStep; } } -mitk::Image::Pointer IgnorePixelMaskGenerator::GetMask() +mitk::Image::ConstPointer IgnorePixelMaskGenerator::GetMask() { if (IsUpdateRequired()) { if (m_inputImage.IsNull()) { MITK_ERROR << "Image not set!"; } if (m_IgnoredPixelValue == std::numeric_limits::min()) { MITK_ERROR << "IgnotePixelValue not set!"; } if (m_TimeStep > (m_inputImage->GetTimeSteps() - 1)) { MITK_ERROR << "Invalid time step: " << m_TimeStep << ". The image has " << m_inputImage->GetTimeSteps() << " timeSteps!"; } // extractimage time slice ImageTimeSelector::Pointer imgTimeSel = ImageTimeSelector::New(); imgTimeSel->SetInput(m_inputImage); imgTimeSel->SetTimeNr(m_TimeStep); imgTimeSel->UpdateLargestPossibleRegion(); mitk::Image::Pointer timeSliceImage = imgTimeSel->GetOutput(); // update m_InternalMask AccessByItk(timeSliceImage, InternalCalculateMask); m_InternalMask->SetGeometry(timeSliceImage->GetGeometry()); this->Modified(); } m_InternalMaskUpdateTime = m_InternalMask->GetMTime(); return m_InternalMask; } template void IgnorePixelMaskGenerator::InternalCalculateMask(typename itk::Image* image) { typedef itk::Image ImageType; typedef itk::Image MaskType; typename MaskType::Pointer mask = MaskType::New(); mask->SetOrigin(image->GetOrigin()); mask->SetSpacing(image->GetSpacing()); mask->SetLargestPossibleRegion(image->GetLargestPossibleRegion()); mask->SetBufferedRegion(image->GetBufferedRegion()); mask->SetDirection(image->GetDirection()); mask->SetNumberOfComponentsPerPixel(image->GetNumberOfComponentsPerPixel()); mask->Allocate(); mask->FillBuffer(1); // iterate over image and mask and set mask=1 if image=m_IgnorePixelValue itk::ImageRegionConstIterator imageIterator(image, image->GetLargestPossibleRegion()); itk::ImageRegionIterator maskIterator(mask, mask->GetLargestPossibleRegion()); for (imageIterator.GoToBegin(); !imageIterator.IsAtEnd(); ++imageIterator, ++maskIterator) { if (imageIterator.Value() == static_cast(m_IgnoredPixelValue)) { maskIterator.Set(0); } } m_InternalMask = GrabItkImageMemory(mask); } bool IgnorePixelMaskGenerator::IsUpdateRequired() const { unsigned long thisClassTimeStamp = this->GetMTime(); unsigned long internalMaskTimeStamp = m_InternalMask->GetMTime(); unsigned long inputImageTimeStamp = m_inputImage->GetMTime(); if (thisClassTimeStamp > m_InternalMaskUpdateTime) // inputs have changed { return true; } if (m_InternalMaskUpdateTime < inputImageTimeStamp) // mask image has changed outside of this class { return true; } if (internalMaskTimeStamp > m_InternalMaskUpdateTime) // internal mask has been changed outside of this class { return true; } return false; } } // end namespace diff --git a/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h b/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h index 08b70337fd..550405c385 100644 --- a/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h +++ b/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h @@ -1,83 +1,84 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef MITKIGNOREPIXELMASKGEN_ #define MITKIGNOREPIXELMASKGEN_ #include #include #include #include #include namespace mitk { /** * @brief The IgnorePixelMaskGenerator class is used to generate a mask that is zero for specific pixel values in the input image. This class requires an input image. */ 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); /** * @brief The mask will be 0 there inputImage==pixelValue and 1 otherwise */ void SetIgnoredPixelValue(RealType pixelValue); /** * @brief Computes and returns the mask */ - mitk::Image::Pointer GetMask() override; + mitk::Image::ConstPointer GetMask() override; /** * @brief SetTimeStep is used to set the time step for which the mask is to be generated * @param timeStep */ void SetTimeStep(unsigned int timeStep) override; protected: IgnorePixelMaskGenerator(): m_IgnoredPixelValue(std::numeric_limits::min()) { m_TimeStep = 0; m_InternalMaskUpdateTime = 0; m_InternalMask = mitk::Image::New(); } ~IgnorePixelMaskGenerator() override{} template void InternalCalculateMask(typename itk::Image* image); private: bool IsUpdateRequired() const; + mitk::Image::Pointer m_InternalMask; RealType m_IgnoredPixelValue; unsigned long m_InternalMaskUpdateTime; }; } #endif diff --git a/Modules/ImageStatistics/mitkImageMaskGenerator.cpp b/Modules/ImageStatistics/mitkImageMaskGenerator.cpp index d8f37a4e21..abf894aca5 100644 --- a/Modules/ImageStatistics/mitkImageMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkImageMaskGenerator.cpp @@ -1,98 +1,99 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include #include #include namespace mitk { -void ImageMaskGenerator::SetImageMask(Image::Pointer maskImage) +void ImageMaskGenerator::SetImageMask(const Image* maskImage) { - if (m_internalMaskImage != maskImage) + if (m_InternalMaskImage != maskImage) { - m_internalMaskImage = maskImage; + m_InternalMaskImage = maskImage; this->Modified(); } } void ImageMaskGenerator::SetTimeStep(unsigned int timeStep) { if (timeStep != m_TimeStep) { m_TimeStep = timeStep; UpdateInternalMask(); } } void ImageMaskGenerator::UpdateInternalMask() { - unsigned int timeStepForExtraction; - - if (m_TimeStep >= m_internalMaskImage->GetTimeSteps()) - { - MITK_WARN << "Warning: time step > number of time steps in mask image, using last time step"; - timeStepForExtraction = m_internalMaskImage->GetTimeSteps() - 1; - } - else - { - timeStepForExtraction = m_TimeStep; - } - ImageTimeSelector::Pointer imageTimeSelector = ImageTimeSelector::New(); - imageTimeSelector->SetInput(m_internalMaskImage); - imageTimeSelector->SetTimeNr(timeStepForExtraction); - imageTimeSelector->UpdateLargestPossibleRegion(); - - m_InternalMask = mitk::Image::New(); - m_InternalMask = imageTimeSelector->GetOutput(); + if (this->m_inputImage.IsNull()) + { + mitkThrow() << "Cannot update internal mask. Input image is not set."; + } + + const auto timeGeo = this->m_inputImage->GetTimeGeometry(); + if (!timeGeo->IsValidTimeStep(this->m_TimeStep)) + { + mitkThrow() << "Cannot update internal mask. Time step selected that is not supported by input image."; + } + + auto timePoint = this->m_inputImage->GetTimeGeometry()->TimeStepToTimePoint(this->m_TimeStep); + m_InternalMask = SelectImageByTimePoint(m_InternalMaskImage, timePoint); + + if (m_InternalMask.IsNull()) + { + MITK_WARN << "Warning: time step > number of time steps in mask image, using last time step"; + m_InternalMask = SelectImageByTimeStep(m_InternalMaskImage, m_InternalMaskImage->GetTimeSteps()-1); + } } -mitk::Image::Pointer ImageMaskGenerator::GetMask() +mitk::Image::ConstPointer ImageMaskGenerator::GetMask() { - if (m_internalMaskImage.IsNull()) + if (m_InternalMaskImage.IsNull()) { MITK_ERROR << "Mask Image is nullptr"; } if (IsUpdateRequired()) { UpdateInternalMask(); } return m_InternalMask; } bool ImageMaskGenerator::IsUpdateRequired() const { unsigned long internalMaskTimeStamp = m_InternalMask->GetMTime(); - unsigned long maskImageTimeStamp = m_internalMaskImage->GetMTime(); + unsigned long maskImageTimeStamp = m_InternalMaskImage->GetMTime(); if (maskImageTimeStamp > internalMaskTimeStamp) // inputs have changed { return true; } if (this->GetMTime() > maskImageTimeStamp) // input has changed { return true; } return false; } } diff --git a/Modules/ImageStatistics/mitkImageMaskGenerator.h b/Modules/ImageStatistics/mitkImageMaskGenerator.h index eed0de6aa6..d42e93da81 100644 --- a/Modules/ImageStatistics/mitkImageMaskGenerator.h +++ b/Modules/ImageStatistics/mitkImageMaskGenerator.h @@ -1,61 +1,62 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkBinaryMaskGenerator #define mitkBinaryMaskGenerator #include #include #include #include #include namespace mitk { class MITKIMAGESTATISTICS_EXPORT ImageMaskGenerator: public MaskGenerator { public: /** Standard Self typedef */ typedef ImageMaskGenerator 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(BinaryImageMaskGenerator, MaskGenerator); - mitk::Image::Pointer GetMask() override; + mitk::Image::ConstPointer GetMask() override; void SetTimeStep(unsigned int timeStep) override; - void SetImageMask(mitk::Image::Pointer maskImage); + void SetImageMask(const mitk::Image* maskImage); protected: ImageMaskGenerator():Superclass(){ m_InternalMaskUpdateTime = 0; m_InternalMask = mitk::Image::New(); } private: bool IsUpdateRequired() const; void UpdateInternalMask(); - mitk::Image::Pointer m_internalMaskImage; + mitk::Image::ConstPointer m_InternalMaskImage; + mitk::Image::ConstPointer m_InternalMask; unsigned long m_InternalMaskUpdateTime; }; } #endif diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index 45720f16a3..6c3f8d916b 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,548 +1,550 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "mitkImageStatisticsCalculator.h" #include #include #include #include #include #include #include #include #include #include #include #include namespace mitk { void ImageStatisticsCalculator::SetInputImage(const mitk::Image *image) { if (image != m_Image) { m_Image = image; this->Modified(); } } void ImageStatisticsCalculator::SetMask(mitk::MaskGenerator *mask) { if (mask != m_MaskGenerator) { m_MaskGenerator = mask; this->Modified(); } } void ImageStatisticsCalculator::SetSecondaryMask(mitk::MaskGenerator *mask) { if (mask != m_SecondaryMaskGenerator) { m_SecondaryMaskGenerator = mask; this->Modified(); } } void ImageStatisticsCalculator::SetNBinsForHistogramStatistics(unsigned int nBins) { if (nBins != m_nBinsForHistogramStatistics) { m_nBinsForHistogramStatistics = nBins; this->Modified(); this->m_UseBinSizeOverNBins = false; } if (m_UseBinSizeOverNBins) { this->Modified(); this->m_UseBinSizeOverNBins = false; } } unsigned int ImageStatisticsCalculator::GetNBinsForHistogramStatistics() const { return m_nBinsForHistogramStatistics; } void ImageStatisticsCalculator::SetBinSizeForHistogramStatistics(double binSize) { if (binSize != m_binSizeForHistogramStatistics) { m_binSizeForHistogramStatistics = binSize; this->Modified(); this->m_UseBinSizeOverNBins = true; } if (!m_UseBinSizeOverNBins) { this->Modified(); this->m_UseBinSizeOverNBins = true; } } double ImageStatisticsCalculator::GetBinSizeForHistogramStatistics() const { return m_binSizeForHistogramStatistics; } mitk::ImageStatisticsContainer* ImageStatisticsCalculator::GetStatistics(LabelIndex label) { if (m_Image.IsNull()) { mitkThrow() << "no image"; } if (!m_Image->IsInitialized()) { mitkThrow() << "Image not initialized!"; } if (IsUpdateRequired(label)) { auto timeGeometry = m_Image->GetTimeGeometry(); // always compute statistics on all timesteps for (unsigned int timeStep = 0; timeStep < m_Image->GetTimeSteps(); timeStep++) { if (m_MaskGenerator.IsNotNull()) { m_MaskGenerator->SetTimeStep(timeStep); //See T25625: otherwise, the mask is not computed again after setting a different time step m_MaskGenerator->Modified(); m_InternalMask = m_MaskGenerator->GetMask(); if (m_MaskGenerator->GetReferenceImage().IsNotNull()) { m_InternalImageForStatistics = m_MaskGenerator->GetReferenceImage(); } else { m_InternalImageForStatistics = m_Image; } } else { m_InternalImageForStatistics = m_Image; } if (m_SecondaryMaskGenerator.IsNotNull()) { m_SecondaryMaskGenerator->SetTimeStep(timeStep); m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(); } ImageTimeSelector::Pointer imgTimeSel = ImageTimeSelector::New(); imgTimeSel->SetInput(m_InternalImageForStatistics); imgTimeSel->SetTimeNr(timeStep); imgTimeSel->UpdateLargestPossibleRegion(); imgTimeSel->Update(); m_ImageTimeSlice = imgTimeSel->GetOutput(); // Calculate statistics with/without mask if (m_MaskGenerator.IsNull() && m_SecondaryMaskGenerator.IsNull()) { // 1) calculate statistics unmasked: AccessByItk_2(m_ImageTimeSlice, InternalCalculateStatisticsUnmasked, timeGeometry, timeStep) } else { // 2) calculate statistics masked AccessByItk_2(m_ImageTimeSlice, InternalCalculateStatisticsMasked, timeGeometry, timeStep) } } } auto it = m_StatisticContainers.find(label); if (it != m_StatisticContainers.end()) { return (it->second).GetPointer(); } else { mitkThrow() << "unknown label"; return nullptr; } } template void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked( typename itk::Image *image, const TimeGeometry *timeGeometry, TimeStepType timeStep) { typedef typename itk::Image ImageType; typedef typename mitk::StatisticsImageFilter ImageStatisticsFilterType; typedef typename itk::MinMaxImageFilterWithIndex MinMaxFilterType; // reset statistics container if exists ImageStatisticsContainer::Pointer statisticContainerForImage; LabelIndex labelNoMask = 1; auto it = m_StatisticContainers.find(labelNoMask); if (it != m_StatisticContainers.end()) { statisticContainerForImage = it->second; } else { statisticContainerForImage = ImageStatisticsContainer::New(); statisticContainerForImage->SetTimeGeometry(const_cast(timeGeometry)); m_StatisticContainers.emplace(labelNoMask, statisticContainerForImage); } auto statObj = ImageStatisticsContainer::ImageStatisticsObject(); typename ImageStatisticsFilterType::Pointer statisticsFilter = ImageStatisticsFilterType::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 MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); minMaxFilter->SetInput(image); minMaxFilter->UpdateLargestPossibleRegion(); typename ImageType::PixelType minval = minMaxFilter->GetMin(); typename ImageType::PixelType maxval = minMaxFilter->GetMax(); typename ImageType::IndexType tmpMinIndex = minMaxFilter->GetMinIndex(); typename ImageType::IndexType tmpMaxIndex = minMaxFilter->GetMaxIndex(); // 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]; } statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUMPOSITION(), minIndex); statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUMPOSITION(), maxIndex); // convert m_binSize in m_nBins if necessary unsigned int nBinsForHistogram; if (m_UseBinSizeOverNBins) { nBinsForHistogram = std::max(static_cast(std::ceil(maxval - minval)) / m_binSizeForHistogramStatistics, 10.); // do not allow less than 10 bins } else { nBinsForHistogram = m_nBinsForHistogramStatistics; } statisticsFilter->SetHistogramParameters(nBinsForHistogram, minval, maxval); try { statisticsFilter->Update(); } catch (const itk::ExceptionObject &e) { mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } auto voxelVolume = GetVoxelVolume(image); auto numberOfPixels = image->GetLargestPossibleRegion().GetNumberOfPixels(); auto volume = static_cast(numberOfPixels) * voxelVolume; auto variance = statisticsFilter->GetSigma() * statisticsFilter->GetSigma(); auto rms = std::sqrt(std::pow(statisticsFilter->GetMean(), 2.) + statisticsFilter->GetVariance()); // variance = sigma^2 statObj.AddStatistic(mitk::ImageStatisticsConstants::NUMBEROFVOXELS(), static_cast(numberOfPixels)); statObj.AddStatistic(mitk::ImageStatisticsConstants::VOLUME(), volume); statObj.AddStatistic(mitk::ImageStatisticsConstants::MEAN(), statisticsFilter->GetMean()); statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUM(), static_cast(statisticsFilter->GetMinimum())); statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUM(), static_cast(statisticsFilter->GetMaximum())); statObj.AddStatistic(mitk::ImageStatisticsConstants::STANDARDDEVIATION(), statisticsFilter->GetSigma()); statObj.AddStatistic(mitk::ImageStatisticsConstants::VARIANCE(), variance); statObj.AddStatistic(mitk::ImageStatisticsConstants::SKEWNESS(), statisticsFilter->GetSkewness()); statObj.AddStatistic(mitk::ImageStatisticsConstants::KURTOSIS(), statisticsFilter->GetKurtosis()); statObj.AddStatistic(mitk::ImageStatisticsConstants::RMS(), rms); statObj.AddStatistic(mitk::ImageStatisticsConstants::MPP(), statisticsFilter->GetMPP()); statObj.AddStatistic(mitk::ImageStatisticsConstants::ENTROPY(), statisticsFilter->GetEntropy()); statObj.AddStatistic(mitk::ImageStatisticsConstants::MEDIAN(), statisticsFilter->GetMedian()); statObj.AddStatistic(mitk::ImageStatisticsConstants::UNIFORMITY(), statisticsFilter->GetUniformity()); statObj.AddStatistic(mitk::ImageStatisticsConstants::UPP(), statisticsFilter->GetUPP()); statObj.m_Histogram = statisticsFilter->GetHistogram(); statisticContainerForImage->SetStatisticsForTimeStep(timeStep, statObj); } template double ImageStatisticsCalculator::GetVoxelVolume(typename itk::Image *image) const { auto spacing = image->GetSpacing(); double voxelVolume = 1.; for (unsigned int i = 0; i < image->GetImageDimension(); i++) { voxelVolume *= spacing[i]; } return voxelVolume; } template void ImageStatisticsCalculator::InternalCalculateStatisticsMasked(typename itk::Image *image, const TimeGeometry *timeGeometry, unsigned int timeStep) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef typename MaskType::PixelType LabelPixelType; typedef LabelStatisticsImageFilter ImageStatisticsFilterType; typedef MaskUtilities MaskUtilType; typedef typename itk::MinMaxLabelImageFilterWithIndex MinMaxLabelFilterType; // workaround: if m_SecondaryMaskGenerator ist not null but m_MaskGenerator is! (this is the case if we request a // 'ignore zuero valued pixels' mask in the gui but do not define a primary mask) bool swapMasks = false; if (m_SecondaryMask.IsNotNull() && m_InternalMask.IsNull()) { m_InternalMask = m_SecondaryMask; m_SecondaryMask = nullptr; swapMasks = true; } // maskImage has to have the same dimension as image - typename MaskType::Pointer maskImage = MaskType::New(); + typename MaskType::ConstPointer maskImage = MaskType::New(); try { // try to access the pixel values directly (no copying or casting). Only works if mask pixels are of pixelType // unsigned short maskImage = ImageToItkImage(m_InternalMask); } catch (const itk::ExceptionObject &) { + typename MaskType::Pointer noneConstMaskImage; //needed to work arround the fact that CastToItkImage currently does not support const itk images. // if the pixel type of the mask is not short, then we have to make a copy of m_InternalMask (and cast the values) - CastToItkImage(m_InternalMask, maskImage); + CastToItkImage(m_InternalMask, noneConstMaskImage); + maskImage = noneConstMaskImage; } // if we have a secondary mask (say a ignoreZeroPixelMask) we need to combine the masks (corresponds to AND) if (m_SecondaryMask.IsNotNull()) { // dirty workaround for a bug when pf mask + any other mask is used in conjunction. We need a proper fix for this // (Fabian Isensee is responsible and probably working on it!) if (m_InternalMask->GetDimension() == 2 && (m_SecondaryMask->GetDimension() == 3 || m_SecondaryMask->GetDimension() == 4)) { mitk::Image::ConstPointer old_img = m_SecondaryMaskGenerator->GetReferenceImage(); m_SecondaryMaskGenerator->SetInputImage(m_MaskGenerator->GetReferenceImage()); m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(); m_SecondaryMaskGenerator->SetInputImage(old_img); } - typename MaskType::Pointer secondaryMaskImage = MaskType::New(); + typename MaskType::ConstPointer secondaryMaskImage = MaskType::New(); secondaryMaskImage = ImageToItkImage(m_SecondaryMask); // secondary mask should be a ignore zero value pixel mask derived from image. it has to be cropped to the mask // region (which may be planar or simply smaller) typename MaskUtilities::Pointer secondaryMaskMaskUtil = MaskUtilities::New(); secondaryMaskMaskUtil->SetImage(secondaryMaskImage.GetPointer()); secondaryMaskMaskUtil->SetMask(maskImage.GetPointer()); - typename MaskType::Pointer adaptedSecondaryMaskImage = secondaryMaskMaskUtil->ExtractMaskImageRegion(); + typename MaskType::ConstPointer adaptedSecondaryMaskImage = secondaryMaskMaskUtil->ExtractMaskImageRegion(); typename itk::MaskImageFilter2::Pointer maskFilter = itk::MaskImageFilter2::New(); maskFilter->SetInput1(maskImage); maskFilter->SetInput2(adaptedSecondaryMaskImage); maskFilter->SetMaskingValue( 1); // all pixels of maskImage where secondaryMaskImage==1 will be kept, all the others are set to 0 maskFilter->UpdateLargestPossibleRegion(); maskImage = maskFilter->GetOutput(); } 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(); + typename ImageType::ConstPointer adaptedImage = ImageType::New(); adaptedImage = maskUtil->ExtractMaskImageRegion(); // this also checks mask sanity // find min, max, minindex and maxindex typename MinMaxLabelFilterType::Pointer minMaxFilter = MinMaxLabelFilterType::New(); minMaxFilter->SetInput(adaptedImage); minMaxFilter->SetLabelInput(maskImage); minMaxFilter->UpdateLargestPossibleRegion(); // set histogram parameters for each label individually (min/max may be different for each label) typedef typename std::unordered_map MapType; std::vector relevantLabels = minMaxFilter->GetRelevantLabels(); MapType minVals; MapType maxVals; std::unordered_map nBins; for (LabelPixelType label : relevantLabels) { minVals[label] = static_cast(minMaxFilter->GetMin(label)); maxVals[label] = static_cast(minMaxFilter->GetMax(label)); unsigned int nBinsForHistogram; if (m_UseBinSizeOverNBins) { nBinsForHistogram = std::max(static_cast(std::ceil(minMaxFilter->GetMax(label) - minMaxFilter->GetMin(label))) / m_binSizeForHistogramStatistics, 10.); // do not allow less than 10 bins } else { nBinsForHistogram = m_nBinsForHistogramStatistics; } nBins[label] = nBinsForHistogram; } typename ImageStatisticsFilterType::Pointer imageStatisticsFilter = ImageStatisticsFilterType::New(); imageStatisticsFilter->SetDirectionTolerance(0.001); imageStatisticsFilter->SetCoordinateTolerance(0.001); imageStatisticsFilter->SetInput(adaptedImage); imageStatisticsFilter->SetLabelInput(maskImage); imageStatisticsFilter->SetHistogramParameters(nBins, minVals, maxVals); imageStatisticsFilter->Update(); auto labels = imageStatisticsFilter->GetValidLabelValues(); auto it = labels.begin(); while (it != labels.end()) { ImageStatisticsContainer::Pointer statisticContainerForLabelImage; auto labelIt = m_StatisticContainers.find(*it); // reset if statisticContainer already exist if (labelIt != m_StatisticContainers.end()) { statisticContainerForLabelImage = labelIt->second; } // create new statisticContainer else { statisticContainerForLabelImage = ImageStatisticsContainer::New(); statisticContainerForLabelImage->SetTimeGeometry(const_cast(timeGeometry)); // link label (*it) to statisticContainer m_StatisticContainers.emplace(*it, statisticContainerForLabelImage); } ImageStatisticsContainer::ImageStatisticsObject statObj; // find min, max, minindex and maxindex // make sure to only look in the masked region, use a masker for this vnl_vector minIndex, maxIndex; mitk::Point3D worldCoordinateMin; mitk::Point3D worldCoordinateMax; mitk::Point3D indexCoordinateMin; mitk::Point3D indexCoordinateMax; m_InternalImageForStatistics->GetGeometry()->IndexToWorld(minMaxFilter->GetMinIndex(*it), worldCoordinateMin); m_InternalImageForStatistics->GetGeometry()->IndexToWorld(minMaxFilter->GetMaxIndex(*it), worldCoordinateMax); m_Image->GetGeometry()->WorldToIndex(worldCoordinateMin, indexCoordinateMin); m_Image->GetGeometry()->WorldToIndex(worldCoordinateMax, indexCoordinateMax); minIndex.set_size(3); maxIndex.set_size(3); // for (unsigned int i=0; i < tmpMaxIndex.GetIndexDimension(); i++) for (unsigned int i = 0; i < 3; i++) { minIndex[i] = indexCoordinateMin[i]; maxIndex[i] = indexCoordinateMax[i]; } statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUMPOSITION(), minIndex); statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUMPOSITION(), maxIndex); auto voxelVolume = GetVoxelVolume(image); auto numberOfVoxels = static_cast(imageStatisticsFilter->GetCount(*it)); auto volume = static_cast(numberOfVoxels) * voxelVolume; auto rms = std::sqrt(std::pow(imageStatisticsFilter->GetMean(*it), 2.) + imageStatisticsFilter->GetVariance(*it)); // variance = sigma^2 auto variance = imageStatisticsFilter->GetSigma(*it) * imageStatisticsFilter->GetSigma(*it); statObj.AddStatistic(mitk::ImageStatisticsConstants::NUMBEROFVOXELS(), numberOfVoxels); statObj.AddStatistic(mitk::ImageStatisticsConstants::VOLUME(), volume); statObj.AddStatistic(mitk::ImageStatisticsConstants::MEAN(), imageStatisticsFilter->GetMean(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUM(), static_cast(imageStatisticsFilter->GetMinimum(*it))); statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUM(), static_cast(imageStatisticsFilter->GetMaximum(*it))); statObj.AddStatistic(mitk::ImageStatisticsConstants::STANDARDDEVIATION(), imageStatisticsFilter->GetSigma(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::VARIANCE(), variance); statObj.AddStatistic(mitk::ImageStatisticsConstants::SKEWNESS(), imageStatisticsFilter->GetSkewness(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::KURTOSIS(), imageStatisticsFilter->GetKurtosis(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::RMS(), rms); statObj.AddStatistic(mitk::ImageStatisticsConstants::MPP(), imageStatisticsFilter->GetMPP(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::ENTROPY(), imageStatisticsFilter->GetEntropy(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::MEDIAN(), imageStatisticsFilter->GetMedian(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::UNIFORMITY(), imageStatisticsFilter->GetUniformity(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::UPP(), imageStatisticsFilter->GetUPP(*it)); statObj.m_Histogram = imageStatisticsFilter->GetHistogram(*it); statisticContainerForLabelImage->SetStatisticsForTimeStep(timeStep, statObj); ++it; } // swap maskGenerators back if (swapMasks) { m_SecondaryMask = m_InternalMask; m_InternalMask = nullptr; } } bool ImageStatisticsCalculator::IsUpdateRequired(LabelIndex label) const { unsigned long thisClassTimeStamp = this->GetMTime(); unsigned long inputImageTimeStamp = m_Image->GetMTime(); auto it = m_StatisticContainers.find(label); if (it == m_StatisticContainers.end()) { return true; } unsigned long statisticsTimeStamp = it->second->GetMTime(); if (thisClassTimeStamp > statisticsTimeStamp) // inputs have changed { return true; } if (inputImageTimeStamp > statisticsTimeStamp) // image has changed { return true; } if (m_MaskGenerator.IsNotNull()) { unsigned long maskGeneratorTimeStamp = m_MaskGenerator->GetMTime(); if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a mask generator and it has changed { return true; } } if (m_SecondaryMaskGenerator.IsNotNull()) { unsigned long maskGeneratorTimeStamp = m_SecondaryMaskGenerator->GetMTime(); if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a secondary mask generator and it has changed { return true; } } return false; } } // namespace mitk diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h index b3489707c5..459aa688e2 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h @@ -1,120 +1,120 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef MITKIMAGESTATISTICSCALCULATOR #define MITKIMAGESTATISTICSCALCULATOR #include #include #include #include namespace mitk { class MITKIMAGESTATISTICS_EXPORT ImageStatisticsCalculator: public itk::Object { public: /** Standard Self typedef */ typedef ImageStatisticsCalculator 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; typedef unsigned short MaskPixelType; using LabelIndex = ImageStatisticsContainer::LabelIndex; /**Documentation @brief Set the image for which the statistics are to be computed.*/ void SetInputImage(const mitk::Image* image); /**Documentation @brief Set the mask generator that creates the mask which is to be used to calculate statistics. If no more mask is desired simply set @param mask to nullptr*/ void SetMask(mitk::MaskGenerator* mask); /**Documentation @brief Set this if more than one mask should be applied (for instance if a IgnorePixelValueMask were to be used alongside with a segmentation). Both masks are combined using pixel wise AND operation. The secondary mask does not have to be the same size than the primary but they need to have some overlap*/ void SetSecondaryMask(mitk::MaskGenerator* mask); /**Documentation @brief Set number of bins to be used for histogram statistics. If Bin size is set after number of bins, bin size will be used instead!*/ void SetNBinsForHistogramStatistics(unsigned int nBins); /**Documentation @brief Retrieve the number of bins used for histogram statistics. Careful: The return value does not indicate whether NBins or BinSize is used. That solely depends on which parameter has been set last.*/ unsigned int GetNBinsForHistogramStatistics() const; /**Documentation @brief Set bin size to be used for histogram statistics. If nbins is set after bin size, nbins will be used instead!*/ void SetBinSizeForHistogramStatistics(double binSize); /**Documentation @brief Retrieve the bin size for histogram statistics. Careful: The return value does not indicate whether NBins or BinSize is used. That solely depends on which parameter has been set last.*/ double GetBinSizeForHistogramStatistics() const; /**Documentation @brief Returns the statistics for label @a label. If these requested statistics are not computed yet the computation is done as well. For performance reasons, statistics for all labels in the image are computed at once. */ ImageStatisticsContainer* GetStatistics(LabelIndex label=1); protected: ImageStatisticsCalculator(){ m_nBinsForHistogramStatistics = 100; m_binSizeForHistogramStatistics = 10; m_UseBinSizeOverNBins = false; }; private: //Calculates statistics for each timestep for image template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsUnmasked( typename itk::Image< TPixel, VImageDimension >* image, const TimeGeometry* timeGeometry, TimeStepType timeStep); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsMasked( typename itk::Image< TPixel, VImageDimension >* image, const TimeGeometry* timeGeometry, unsigned int timeStep); template < typename TPixel, unsigned int VImageDimension > double GetVoxelVolume(typename itk::Image* image) const; bool IsUpdateRequired(LabelIndex label) const; mitk::Image::ConstPointer m_Image; mitk::Image::Pointer m_ImageTimeSlice; mitk::Image::ConstPointer m_InternalImageForStatistics; mitk::MaskGenerator::Pointer m_MaskGenerator; - mitk::Image::Pointer m_InternalMask; + mitk::Image::ConstPointer m_InternalMask; mitk::MaskGenerator::Pointer m_SecondaryMaskGenerator; - mitk::Image::Pointer m_SecondaryMask; + mitk::Image::ConstPointer m_SecondaryMask; unsigned int m_nBinsForHistogramStatistics; double m_binSizeForHistogramStatistics; bool m_UseBinSizeOverNBins; std::map m_StatisticContainers; }; } #endif // MITKIMAGESTATISTICSCALCULATOR diff --git a/Modules/ImageStatistics/mitkMaskGenerator.cpp b/Modules/ImageStatistics/mitkMaskGenerator.cpp index af8b4a812b..0e29bbb2e4 100644 --- a/Modules/ImageStatistics/mitkMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkMaskGenerator.cpp @@ -1,64 +1,45 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include namespace mitk { MaskGenerator::MaskGenerator(): m_TimeStep(0) { m_inputImage = nullptr; } -mitk::Image::Pointer MaskGenerator::GetMask() -{ - return mitk::Image::New(); -} - - -//typename itk::Region<3>::Pointer MaskGenerator::GetImageRegionOfMask(Image::Pointer image) -//{ -// if (m_InternalMask.IsNull() || m_Modified) -// { -// MITK_ERROR << "Update MaskGenerator first!"; -// } - -// mitk::BaseGeometry::Pointer imageGeometry = image->GetGeometry(); -// mitk::BaseGeometry::Pointer maskGeometry = m_InternalMask->GetGeometry(); - - -//} - void MaskGenerator::SetTimeStep(unsigned int timeStep) { if (timeStep != m_TimeStep) { m_TimeStep = timeStep; } } void MaskGenerator::SetInputImage(mitk::Image::ConstPointer inputImg) { if (inputImg != m_inputImage) { m_inputImage = inputImg; this->Modified(); } } mitk::Image::ConstPointer MaskGenerator::GetReferenceImage() { return m_inputImage; } } diff --git a/Modules/ImageStatistics/mitkMaskGenerator.h b/Modules/ImageStatistics/mitkMaskGenerator.h index 76d976c777..adbafda62c 100644 --- a/Modules/ImageStatistics/mitkMaskGenerator.h +++ b/Modules/ImageStatistics/mitkMaskGenerator.h @@ -1,76 +1,74 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef 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(); + virtual mitk::Image::ConstPointer GetMask() = 0; /** * @brief GetReferenceImage per default returns the inputImage (as set by SetInputImage). If no input image is set it will return a nullptr. */ virtual mitk::Image::ConstPointer GetReferenceImage(); /** * @brief SetInputImage is used to set the input image to the mask generator. Some subclasses require an input image, others don't. See the documentation of the specific Mask Generator for more information. */ void SetInputImage(mitk::Image::ConstPointer inputImg); virtual void SetTimeStep(unsigned int timeStep); protected: MaskGenerator(); unsigned int m_TimeStep; - mitk::Image::Pointer m_InternalMask; mitk::Image::ConstPointer m_inputImage; private: }; } #endif // MITKMASKGENERATOR diff --git a/Modules/ImageStatistics/mitkMaskUtilities.h b/Modules/ImageStatistics/mitkMaskUtilities.h index ec0ff30f99..4ec2b2631e 100644 --- a/Modules/ImageStatistics/mitkMaskUtilities.h +++ b/Modules/ImageStatistics/mitkMaskUtilities.h @@ -1,88 +1,88 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef MITKMASKUTIL #define MITKMASKUTIL #include #include #include #include namespace mitk { /** * @brief Utility class for mask operations. It checks whether an image and a mask are compatible (spacing, orientation, etc...) * and it can also crop an image to the LargestPossibleRegion of the Mask */ template class MaskUtilities: public itk::Object { public: /** Standard Self typedef */ typedef MaskUtilities 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(MaskUtilities, itk::Object); typedef itk::Image ImageType; typedef itk::Image MaskType; /** * @brief Set image */ - void SetImage(ImageType* image); + void SetImage(const ImageType* image); /** * @brief Set mask */ - void SetMask(MaskType* mask); + void SetMask(const MaskType* mask); /** * @brief Checks whether mask and image are compatible for joint access (as via iterators). * Spacing and direction must be the same between the two and they must be aligned. Also, the mask must be completely inside the image */ bool CheckMaskSanity(); /** * @brief Crops the image to the LargestPossibleRegion of the mask */ - typename itk::Image::Pointer ExtractMaskImageRegion(); + typename ImageType::ConstPointer ExtractMaskImageRegion(); protected: MaskUtilities(): m_Image(nullptr), m_Mask(nullptr){} ~MaskUtilities() override{} private: - itk::Image* m_Image; - itk::Image* m_Mask; + const ImageType* m_Image; + const MaskType* m_Mask; }; /** Tolerance used to check if the mask and input image are compatible for * coordinate aspects (orgin, size, grid alignment).*/ constexpr double MASK_SUITABILITY_TOLERANCE_COORDINATE = NODE_PREDICATE_GEOMETRY_DEFAULT_CHECK_COORDINATE_PRECISION; /** Tolerance used to check if the mask and input image are compatible for * direction aspects (orientation of mask and image).*/ constexpr double MASK_SUITABILITY_TOLERANCE_DIRECTION = NODE_PREDICATE_GEOMETRY_DEFAULT_CHECK_DIRECTION_PRECISION; } #ifndef ITK_MANUAL_INSTANTIATION #include #endif #endif diff --git a/Modules/ImageStatistics/mitkMaskUtilities.tpp b/Modules/ImageStatistics/mitkMaskUtilities.tpp index 25a8b4a2ff..f1d52ea2aa 100644 --- a/Modules/ImageStatistics/mitkMaskUtilities.tpp +++ b/Modules/ImageStatistics/mitkMaskUtilities.tpp @@ -1,196 +1,194 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef MITKMASKUTIL_TPP #define MITKMASKUTIL_TPP #include #include #include #include #include namespace mitk { template - void MaskUtilities::SetImage(ImageType* image) + void MaskUtilities::SetImage(const ImageType* image) { if (image != m_Image) { m_Image = image; } } template - void MaskUtilities::SetMask(MaskType* mask) + void MaskUtilities::SetMask(const MaskType* mask) { if (mask != m_Mask) { m_Mask = mask; } } template bool MaskUtilities::CheckMaskSanity() { if (m_Mask==nullptr || m_Image==nullptr) { MITK_ERROR << "Set an image and a mask first"; } typedef itk::Image< TPixel, VImageDimension > ImageType; typedef typename ImageType::PointType PointType; typedef typename ImageType::DirectionType DirectionType; bool maskSanity = true; if (m_Mask==nullptr) { MITK_ERROR << "Something went wrong when casting the mitk mask image to an itk mask image. Do the mask and the input image have the same dimension?"; // note to self: We could try to convert say a 2d mask to a 3d mask if the image is 3d. (mask and image dimension have to match.) } // check direction DirectionType imageDirection = m_Image->GetDirection(); DirectionType maskDirection = m_Mask->GetDirection(); for(unsigned int i = 0; i < imageDirection.ColumnDimensions; ++i ) { for(unsigned int j = 0; j < imageDirection.ColumnDimensions; ++j ) { double differenceDirection = imageDirection[i][j] - maskDirection[i][j]; if (fabs(differenceDirection) > MASK_SUITABILITY_TOLERANCE_DIRECTION) { maskSanity = false; MITK_INFO << "Mask needs to have same direction as image! (Image direction: " << imageDirection << "; Mask direction: " << maskDirection << ")"; } } } // check spacing PointType imageSpacing = m_Image->GetSpacing(); PointType maskSpacing = m_Mask->GetSpacing(); for (unsigned int i = 0; i < VImageDimension; i++) { if ( fabs( maskSpacing[i] - imageSpacing[i] ) > MASK_SUITABILITY_TOLERANCE_COORDINATE ) { maskSanity = false; MITK_INFO << "Spacing of mask and image is not equal. Mask: " << maskSpacing << " image: " << imageSpacing; } } // check alignment // Make sure that the voxels of mask and image are correctly "aligned", i.e., voxel boundaries are the same in both images PointType imageOrigin = m_Image->GetOrigin(); PointType maskOrigin = m_Mask->GetOrigin(); typedef itk::ContinuousIndex ContinousIndexType; ContinousIndexType maskOriginContinousIndex, imageOriginContinousIndex; m_Image->TransformPhysicalPointToContinuousIndex(maskOrigin, maskOriginContinousIndex); m_Image->TransformPhysicalPointToContinuousIndex(imageOrigin, imageOriginContinousIndex); for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { double misalignment = maskOriginContinousIndex[i] - floor( maskOriginContinousIndex[i] + 0.5 ); // misalignment must be a multiple (int) of spacing in that direction if ( fmod(misalignment,imageSpacing[i]) > MASK_SUITABILITY_TOLERANCE_COORDINATE) { maskSanity = false; MITK_INFO << "Pixels/voxels of mask and image are not sufficiently aligned! (Misalignment: " << fmod(misalignment,imageSpacing[i]) << ")"; } } // mask must be completely inside image region // Make sure that mask region is contained within image region if ( m_Mask!=nullptr && !m_Image->GetLargestPossibleRegion().IsInside( m_Mask->GetLargestPossibleRegion() ) ) { maskSanity = false; MITK_INFO << "Mask region needs to be inside of image region! (Image region: " << m_Image->GetLargestPossibleRegion() << "; Mask region: " << m_Mask->GetLargestPossibleRegion() << ")"; } return maskSanity; } template - typename itk::Image::Pointer MaskUtilities::ExtractMaskImageRegion() + typename MaskUtilities::ImageType::ConstPointer MaskUtilities::ExtractMaskImageRegion() { if (m_Mask==nullptr || m_Image==nullptr) { MITK_ERROR << "Set an image and a mask first"; } bool maskSanity = CheckMaskSanity(); if (!maskSanity) { MITK_ERROR << "Mask and image are not compatible"; } - typedef itk::Image< TPixel, VImageDimension > ImageType; - typedef itk::Image< unsigned short, VImageDimension > MaskType; typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractImageFilterType; typename ImageType::SizeType imageSize = m_Image->GetBufferedRegion().GetSize(); typename ImageType::SizeType maskSize = m_Mask->GetBufferedRegion().GetSize(); - typename itk::Image::Pointer extractedImg = itk::Image::New(); + typename itk::Image::ConstPointer resultImg; bool maskSmallerImage = false; for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { if ( maskSize[i] < imageSize[i] ) { maskSmallerImage = true; } } if ( maskSmallerImage ) { typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New(); typename MaskType::PointType maskOrigin = m_Mask->GetOrigin(); typename ImageType::PointType imageOrigin = m_Image->GetOrigin(); typename MaskType::SpacingType maskSpacing = m_Mask->GetSpacing(); typename ImageType::RegionType extractionRegion; typename ImageType::IndexType extractionRegionIndex; for (unsigned int i=0; i < maskOrigin.GetPointDimension(); i++) { extractionRegionIndex[i] = (maskOrigin[i] - imageOrigin[i]) / maskSpacing[i]; } extractionRegion.SetIndex(extractionRegionIndex); extractionRegion.SetSize(m_Mask->GetLargestPossibleRegion().GetSize()); extractImageFilter->SetInput( m_Image ); extractImageFilter->SetExtractionRegion( extractionRegion ); extractImageFilter->SetCoordinateTolerance(MASK_SUITABILITY_TOLERANCE_COORDINATE); extractImageFilter->SetDirectionTolerance(MASK_SUITABILITY_TOLERANCE_DIRECTION); extractImageFilter->Update(); - extractedImg = extractImageFilter->GetOutput(); + auto extractedImg = extractImageFilter->GetOutput(); extractedImg->SetOrigin(m_Mask->GetOrigin()); extractedImg->SetLargestPossibleRegion(m_Mask->GetLargestPossibleRegion()); extractedImg->SetBufferedRegion(m_Mask->GetBufferedRegion()); - + resultImg = extractedImg; } else { - extractedImg = m_Image; + resultImg = m_Image; } - return extractedImg; + return resultImg; } } #endif diff --git a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp index 6acffbfdf4..c44d2100d6 100644 --- a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp @@ -1,512 +1,512 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include #include #include #include "mitkImageAccessByItk.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include namespace mitk { -void PlanarFigureMaskGenerator::SetPlanarFigure(mitk::PlanarFigure::Pointer planarFigure) +void PlanarFigureMaskGenerator::SetPlanarFigure(mitk::PlanarFigure* planarFigure) { - if ( planarFigure.IsNull() ) + if (nullptr == planarFigure ) { throw std::runtime_error( "Error: planar figure empty!" ); } const PlaneGeometry *planarFigurePlaneGeometry = planarFigure->GetPlaneGeometry(); if ( planarFigurePlaneGeometry == nullptr ) { throw std::runtime_error( "Planar-Figure not yet initialized!" ); } const auto *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry ); if ( planarFigureGeometry == nullptr ) { throw std::runtime_error( "Non-planar planar figures not supported!" ); } if (planarFigure != m_PlanarFigure) { this->Modified(); m_PlanarFigure = planarFigure; } } mitk::Image::ConstPointer PlanarFigureMaskGenerator::GetReferenceImage() { if (IsUpdateRequired()) { this->CalculateMask(); } return m_ReferenceImage; } template < typename TPixel, unsigned int VImageDimension > void PlanarFigureMaskGenerator::InternalCalculateMaskFromPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ) { typedef itk::Image< unsigned short, 2 > MaskImage2DType; 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_inputImage->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 vtkSmartPointer points = vtkSmartPointer::New(); for (const auto& point : planarFigurePolyline) { Point3D point3D; // Convert 2D point back to the local index coordinates of the selected image planarFigurePlaneGeometry->Map(point, point3D); imageGeometry3D->WorldToIndex(point3D, point3D); points->InsertNextPoint(point3D[i0], point3D[i1], 0); } vtkSmartPointer holePoints; if (!planarFigureHolePolyline.empty()) { holePoints = vtkSmartPointer::New(); Point3D point3D; for (const auto& point : planarFigureHolePolyline) { planarFigurePlaneGeometry->Map(point, 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}; 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."; } // 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() ); 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(); } template < typename TPixel, unsigned int VImageDimension > void PlanarFigureMaskGenerator::InternalCalculateMaskFromOpenPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ) { typedef itk::Image< unsigned short, 2 > MaskImage2DType; typedef itk::LineIterator< MaskImage2DType > LineIteratorType; typedef MaskImage2DType::IndexType IndexType2D; typedef std::vector< IndexType2D > IndexVecType; 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(0); // all PolylinePoints of the PlanarFigure are stored in a vtkPoints object. const mitk::PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); const mitk::BaseGeometry *imageGeometry3D = m_inputImage->GetGeometry( 0 ); // 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; } int numPolyLines = m_PlanarFigure->GetPolyLinesSize(); for ( int lineId = 0; lineId < numPolyLines; ++lineId ) { // store the polyline contour as vtkPoints object IndexVecType pointIndices; for(const auto& point : planarFigurePolyline) { Point3D point3D; planarFigurePlaneGeometry->Map(point, point3D); imageGeometry3D->WorldToIndex(point3D, point3D); IndexType2D index2D; index2D[0] = point3D[i0]; index2D[1] = point3D[i1]; pointIndices.push_back( index2D ); } size_t numLineSegments = pointIndices.size() - 1; for (size_t i = 0; i < numLineSegments; ++i) { LineIteratorType lineIt(maskImage, pointIndices[i], pointIndices[i+1]); while (!lineIt.IsAtEnd()) { lineIt.Set(1); ++lineIt; } } } // Store mask m_InternalITKImageMask2D = maskImage; } bool PlanarFigureMaskGenerator::CheckPlanarFigureIsNotTilted(const PlaneGeometry* planarGeometry, const BaseGeometry *geometry) { if (!planarGeometry) return false; if (!geometry) return false; unsigned int axis; return GetPrincipalAxis(geometry,planarGeometry->GetNormal(), axis); } 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(); //normal mitk::eps is to pedantic for this check. See e.g. T27122 //therefore choose a larger epsilon. The value was set a) as small as //possible but b) still allowing to datasets like in (T27122) to pass //when floating rounding errors sum up. const double epsilon = 5e-5; if ( fabs( fabs( axisVector * vector ) - 1.0) < epsilon) { axis = i; return true; } } return false; } void PlanarFigureMaskGenerator::CalculateMask() { if (m_inputImage.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)."; } const BaseGeometry *imageGeometry = m_inputImage->GetGeometry(); if ( imageGeometry == nullptr ) { throw std::runtime_error( "Image geometry invalid!" ); } if (m_inputImage->GetTimeSteps() > 0) { mitk::ImageTimeSelector::Pointer imgTimeSel = mitk::ImageTimeSelector::New(); imgTimeSel->SetInput(m_inputImage); imgTimeSel->SetTimeNr(m_TimeStep); imgTimeSel->UpdateLargestPossibleRegion(); m_InternalTimeSliceImage = imgTimeSel->GetOutput(); } else { m_InternalTimeSliceImage = m_inputImage; } m_InternalITKImageMask2D = nullptr; const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); const auto *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry ); //const BaseGeometry *imageGeometry = m_inputImage->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 itk::Image< unsigned short, 3 >::IndexType index; imageGeometry->WorldToIndex( planarFigureGeometry->GetOrigin(), index ); unsigned int slice = index[axis]; m_PlanarFigureSlice = slice; // extract image slice which corresponds to the planarFigure and store it in m_InternalImageSlice mitk::Image::ConstPointer inputImageSlice = extract2DImageSlice(axis, slice); //mitk::IOUtil::Save(inputImageSlice, "/home/fabian/inputSliceImage.nrrd"); // Compute mask from PlanarFigure // rastering for open planar figure: if ( !m_PlanarFigure->IsClosed() ) { AccessFixedDimensionByItk_1(inputImageSlice, InternalCalculateMaskFromOpenPlanarFigure, 2, axis) } else//for closed planar figure { 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); //mitk::IOUtil::Save(planarFigureMaskImage, "/home/fabian/planarFigureMaskImage.nrrd"); //Convert2Dto3DImageFilter::Pointer sliceTo3DImageConverter = Convert2Dto3DImageFilter::New(); //sliceTo3DImageConverter->SetInput(planarFigureMaskImage); //sliceTo3DImageConverter->Update(); //mitk::IOUtil::Save(sliceTo3DImageConverter->GetOutput(), "/home/fabian/3DsliceImage.nrrd"); m_ReferenceImage = inputImageSlice; //mitk::IOUtil::Save(m_ReferenceImage, "/home/fabian/referenceImage.nrrd"); m_InternalMask = planarFigureMaskImage; } void PlanarFigureMaskGenerator::SetTimeStep(unsigned int timeStep) { if (timeStep != m_TimeStep) { m_TimeStep = timeStep; } } -mitk::Image::Pointer PlanarFigureMaskGenerator::GetMask() +mitk::Image::ConstPointer PlanarFigureMaskGenerator::GetMask() { if (IsUpdateRequired()) { this->CalculateMask(); this->Modified(); } m_InternalMaskUpdateTime = this->GetMTime(); return m_InternalMask; } mitk::Image::ConstPointer PlanarFigureMaskGenerator::extract2DImageSlice(unsigned int axis, unsigned int slice) { // Extract slice with given position and direction from image unsigned int dimension = m_InternalTimeSliceImage->GetDimension(); if (dimension == 3) { ExtractImageFilter::Pointer imageExtractor = ExtractImageFilter::New(); imageExtractor->SetInput( m_InternalTimeSliceImage ); imageExtractor->SetSliceDimension( axis ); imageExtractor->SetSliceIndex( slice ); imageExtractor->Update(); return imageExtractor->GetOutput(); } else if(dimension == 2) { return m_InternalTimeSliceImage; } else { MITK_ERROR << "Unsupported image dimension. Dimension is: " << dimension << ". Only 2D and 3D images are supported."; return nullptr; } } bool PlanarFigureMaskGenerator::IsUpdateRequired() const { unsigned long thisClassTimeStamp = this->GetMTime(); unsigned long internalMaskTimeStamp = m_InternalMask->GetMTime(); unsigned long planarFigureTimeStamp = m_PlanarFigure->GetMTime(); unsigned long inputImageTimeStamp = m_inputImage->GetMTime(); if (thisClassTimeStamp > m_InternalMaskUpdateTime) // inputs have changed { return true; } if (m_InternalMaskUpdateTime < planarFigureTimeStamp || m_InternalMaskUpdateTime < inputImageTimeStamp) // mask image has changed outside of this class { return true; } if (internalMaskTimeStamp > m_InternalMaskUpdateTime) // internal mask has been changed outside of this class { return true; } return false; } } diff --git a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h index 0e61304fef..5c2559fc9d 100644 --- a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h +++ b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h @@ -1,145 +1,146 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef MITKPLANARFIGUREMASKGENERATOR #define MITKPLANARFIGUREMASKGENERATOR #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 Pointer; typedef itk::SmartPointer ConstPointer; /** Method for creation through the object factory. */ itkNewMacro(Self); /** Runtime information support. */ itkTypeMacro(PlanarFigureMaskGenerator, MaskGenerator); /** * @brief GetMask Computes and returns the mask * @return mitk::Image::Pointer of the generated mask */ - mitk::Image::Pointer GetMask() override; + mitk::Image::ConstPointer GetMask() override; - void SetPlanarFigure(mitk::PlanarFigure::Pointer planarFigure); + void SetPlanarFigure(mitk::PlanarFigure* planarFigure); mitk::Image::ConstPointer GetReferenceImage() override; /** * @brief SetTimeStep is used to set the time step for which the mask is to be generated * @param timeStep */ void SetTimeStep(unsigned int timeStep) override; itkGetConstMacro(PlanarFigureAxis, unsigned int); itkGetConstMacro(PlanarFigureSlice, unsigned int); /** Helper function that indicates if a passed planar geometry is tilted regarding a given geometry and its main axis. *@pre If either planarGeometry or geometry is nullptr it will return false.*/ static bool CheckPlanarFigureIsNotTilted(const PlaneGeometry* planarGeometry, const BaseGeometry *geometry); protected: PlanarFigureMaskGenerator() : Superclass(), m_ReferenceImage(nullptr), m_PlanarFigureAxis(0), m_InternalMaskUpdateTime(0), m_PlanarFigureSlice(0) { m_InternalMask = mitk::Image::New(); } private: void CalculateMask(); template void InternalCalculateMaskFromPlanarFigure(const itk::Image *image, unsigned int axis); template void InternalCalculateMaskFromOpenPlanarFigure(const itk::Image *image, unsigned int axis); mitk::Image::ConstPointer extract2DImageSlice(unsigned int axis, unsigned int slice); /** Helper function that deduces if the passed vector is equal to one of the primary axis of the geometry.*/ static 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()); } bool IsUpdateRequired() const; mitk::PlanarFigure::Pointer m_PlanarFigure; itk::Image::Pointer m_InternalITKImageMask2D; mitk::Image::ConstPointer m_InternalTimeSliceImage; mitk::Image::ConstPointer m_ReferenceImage; unsigned int m_PlanarFigureAxis; unsigned long m_InternalMaskUpdateTime; unsigned int m_PlanarFigureSlice; + mitk::Image::Pointer m_InternalMask; }; } // namespace mitk #endif // MITKPLANARFIGUREMASKGENERATOR diff --git a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationJob.cpp b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationJob.cpp index c07bcda386..f3cf23d68b 100644 --- a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationJob.cpp +++ b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationJob.cpp @@ -1,221 +1,222 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "QmitkImageStatisticsCalculationJob.h" #include "mitkImageStatisticsCalculator.h" #include #include #include QmitkImageStatisticsCalculationJob::QmitkImageStatisticsCalculationJob() : QThread() , m_StatisticsImage(nullptr) , m_BinaryMask(nullptr) , m_PlanarFigureMask(nullptr) , m_IgnoreZeros(false) , m_HistogramNBins(100) , m_CalculationSuccessful(false) { } QmitkImageStatisticsCalculationJob::~QmitkImageStatisticsCalculationJob() { } void QmitkImageStatisticsCalculationJob::Initialize(const mitk::Image *image, const mitk::Image *binaryImage, const mitk::PlanarFigure *planarFig) { this->m_StatisticsImage = image; this->m_BinaryMask = binaryImage; this->m_PlanarFigureMask = planarFig; } mitk::ImageStatisticsContainer* QmitkImageStatisticsCalculationJob::GetStatisticsData() const { return this->m_StatisticsContainer.GetPointer(); } const mitk::Image* QmitkImageStatisticsCalculationJob::GetStatisticsImage() const { return this->m_StatisticsImage.GetPointer(); } const mitk::Image* QmitkImageStatisticsCalculationJob::GetMaskImage() const { return this->m_BinaryMask.GetPointer(); } const mitk::PlanarFigure* QmitkImageStatisticsCalculationJob::GetPlanarFigure() const { return this->m_PlanarFigureMask.GetPointer(); } void QmitkImageStatisticsCalculationJob::SetIgnoreZeroValueVoxel(bool _arg) { this->m_IgnoreZeros = _arg; } bool QmitkImageStatisticsCalculationJob::GetIgnoreZeroValueVoxel() const { return this->m_IgnoreZeros; } void QmitkImageStatisticsCalculationJob::SetHistogramNBins(unsigned int nbins) { this->m_HistogramNBins = nbins; } unsigned int QmitkImageStatisticsCalculationJob::GetHistogramNBins() const { return this->m_HistogramNBins; } std::string QmitkImageStatisticsCalculationJob::GetLastErrorMessage() const { return m_message; } const QmitkImageStatisticsCalculationJob::HistogramType* QmitkImageStatisticsCalculationJob::GetTimeStepHistogram(unsigned int t) const { if (t >= this->m_HistogramVector.size()) return nullptr; return this->m_HistogramVector.at(t).GetPointer(); } bool QmitkImageStatisticsCalculationJob::GetStatisticsUpdateSuccessFlag() const { return m_CalculationSuccessful; } void QmitkImageStatisticsCalculationJob::run() { bool statisticCalculationSuccessful = true; mitk::ImageStatisticsCalculator::Pointer calculator = mitk::ImageStatisticsCalculator::New(); if(this->m_StatisticsImage.IsNotNull()) { calculator->SetInputImage(m_StatisticsImage); } else { statisticCalculationSuccessful = false; } // Bug 13416 : The ImageStatistics::SetImageMask() method can throw exceptions, i.e. when the dimensionality // of the masked and input image differ, we need to catch them and mark the calculation as failed // the same holds for the ::SetPlanarFigure() try { if(this->m_BinaryMask.IsNotNull()) { mitk::ImageMaskGenerator::Pointer imgMask = mitk::ImageMaskGenerator::New(); - imgMask->SetImageMask(m_BinaryMask->Clone()); + imgMask->SetInputImage(m_StatisticsImage); + imgMask->SetImageMask(m_BinaryMask); calculator->SetMask(imgMask.GetPointer()); } if(this->m_PlanarFigureMask.IsNotNull()) { mitk::PlanarFigureMaskGenerator::Pointer pfMaskGen = mitk::PlanarFigureMaskGenerator::New(); pfMaskGen->SetInputImage(m_StatisticsImage); pfMaskGen->SetPlanarFigure(m_PlanarFigureMask->Clone()); calculator->SetMask(pfMaskGen.GetPointer()); } } catch (const mitk::Exception& e) { MITK_ERROR << "MITK Exception: " << e.what(); m_message = e.what(); statisticCalculationSuccessful = false; } catch( const itk::ExceptionObject& e) { MITK_ERROR << "ITK Exception:" << e.what(); m_message = e.what(); statisticCalculationSuccessful = false; } catch ( const std::runtime_error &e ) { MITK_ERROR<< "Runtime Exception: " << e.what(); m_message = e.what(); statisticCalculationSuccessful = false; } catch ( const std::exception &e ) { MITK_ERROR<< "Standard Exception: " << e.what(); m_message = e.what(); statisticCalculationSuccessful = false; } if (this->m_IgnoreZeros) { mitk::IgnorePixelMaskGenerator::Pointer ignorePixelValueMaskGen = mitk::IgnorePixelMaskGenerator::New(); ignorePixelValueMaskGen->SetIgnoredPixelValue(0); ignorePixelValueMaskGen->SetInputImage(m_StatisticsImage); calculator->SetSecondaryMask(ignorePixelValueMaskGen.GetPointer()); } else { calculator->SetSecondaryMask(nullptr); } calculator->SetNBinsForHistogramStatistics(m_HistogramNBins); try { calculator->GetStatistics(); } catch ( mitk::Exception& e) { m_message = e.GetDescription(); MITK_ERROR<< "MITK Exception: " << e.what(); statisticCalculationSuccessful = false; } catch ( const std::runtime_error &e ) { m_message = "Failure: " + std::string(e.what()); MITK_ERROR<< "Runtime Exception: " << e.what(); statisticCalculationSuccessful = false; } catch ( const std::exception &e ) { m_message = "Failure: " + std::string(e.what()); MITK_ERROR<< "Standard Exception: " << e.what(); statisticCalculationSuccessful = false; } this->m_CalculationSuccessful = statisticCalculationSuccessful; if(statisticCalculationSuccessful) { m_StatisticsContainer = calculator->GetStatistics(); this->m_HistogramVector.clear(); for (unsigned int i = 0; i < m_StatisticsImage->GetTimeSteps(); i++) { HistogramType::ConstPointer tempHistogram; try { if(calculator->GetStatistics()->TimeStepExists(i)) { tempHistogram = calculator->GetStatistics()->GetStatisticsForTimeStep(i).m_Histogram; this->m_HistogramVector.push_back(tempHistogram); } } catch (mitk::Exception&) { MITK_WARN << ":-("; } } } } diff --git a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationRunnable.cpp b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationRunnable.cpp index 51850b06ed..7dc62e5227 100644 --- a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationRunnable.cpp +++ b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationRunnable.cpp @@ -1,186 +1,187 @@ /*=================================================================== 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 "QmitkImageStatisticsCalculationRunnable.h" #include "mitkImageStatisticsCalculator.h" #include #include #include #include "mitkStatisticsToImageRelationRule.h" #include "mitkStatisticsToMaskRelationRule.h" #include "mitkImageStatisticsContainerManager.h" #include "mitkProperties.h" QmitkImageStatisticsCalculationRunnable::QmitkImageStatisticsCalculationRunnable() : QmitkDataGenerationJobBase() , m_StatisticsImage(nullptr) , m_BinaryMask(nullptr) , m_PlanarFigureMask(nullptr) , m_IgnoreZeros(false) , m_HistogramNBins(100) { } QmitkImageStatisticsCalculationRunnable::~QmitkImageStatisticsCalculationRunnable() { } void QmitkImageStatisticsCalculationRunnable::Initialize(const mitk::Image *image, const mitk::Image *binaryImage, const mitk::PlanarFigure *planarFig) { this->m_StatisticsImage = image; this->m_BinaryMask = binaryImage; this->m_PlanarFigureMask = planarFig; } mitk::ImageStatisticsContainer* QmitkImageStatisticsCalculationRunnable::GetStatisticsData() const { return this->m_StatisticsContainer.GetPointer(); } const mitk::Image* QmitkImageStatisticsCalculationRunnable::GetStatisticsImage() const { return this->m_StatisticsImage.GetPointer(); } const mitk::Image* QmitkImageStatisticsCalculationRunnable::GetMaskImage() const { return this->m_BinaryMask.GetPointer(); } const mitk::PlanarFigure* QmitkImageStatisticsCalculationRunnable::GetPlanarFigure() const { return this->m_PlanarFigureMask.GetPointer(); } void QmitkImageStatisticsCalculationRunnable::SetIgnoreZeroValueVoxel(bool _arg) { this->m_IgnoreZeros = _arg; } bool QmitkImageStatisticsCalculationRunnable::GetIgnoreZeroValueVoxel() const { return this->m_IgnoreZeros; } void QmitkImageStatisticsCalculationRunnable::SetHistogramNBins(unsigned int nbins) { this->m_HistogramNBins = nbins; } unsigned int QmitkImageStatisticsCalculationRunnable::GetHistogramNBins() const { return this->m_HistogramNBins; } QmitkDataGenerationJobBase::ResultMapType QmitkImageStatisticsCalculationRunnable::GetResults() const { ResultMapType result; result.emplace("statistics", this->GetStatisticsData()); return result; } bool QmitkImageStatisticsCalculationRunnable::RunComputation() { bool statisticCalculationSuccessful = true; mitk::ImageStatisticsCalculator::Pointer calculator = mitk::ImageStatisticsCalculator::New(); if (this->m_StatisticsImage.IsNotNull()) { calculator->SetInputImage(m_StatisticsImage); } else { statisticCalculationSuccessful = false; } // Bug 13416 : The ImageStatistics::SetImageMask() method can throw exceptions, i.e. when the dimensionality // of the masked and input image differ, we need to catch them and mark the calculation as failed // the same holds for the ::SetPlanarFigure() try { if (this->m_BinaryMask.IsNotNull()) { mitk::ImageMaskGenerator::Pointer imgMask = mitk::ImageMaskGenerator::New(); - imgMask->SetImageMask(m_BinaryMask->Clone()); + imgMask->SetInputImage(m_StatisticsImage); + imgMask->SetImageMask(m_BinaryMask); calculator->SetMask(imgMask.GetPointer()); } if (this->m_PlanarFigureMask.IsNotNull()) { mitk::PlanarFigureMaskGenerator::Pointer pfMaskGen = mitk::PlanarFigureMaskGenerator::New(); pfMaskGen->SetInputImage(m_StatisticsImage); pfMaskGen->SetPlanarFigure(m_PlanarFigureMask->Clone()); calculator->SetMask(pfMaskGen.GetPointer()); } } catch (const std::exception &e) { MITK_ERROR << "Error while configuring the statistics calculator: " << e.what(); m_LastErrorMessage = e.what(); statisticCalculationSuccessful = false; } if (this->m_IgnoreZeros) { mitk::IgnorePixelMaskGenerator::Pointer ignorePixelValueMaskGen = mitk::IgnorePixelMaskGenerator::New(); ignorePixelValueMaskGen->SetIgnoredPixelValue(0); ignorePixelValueMaskGen->SetInputImage(m_StatisticsImage); calculator->SetSecondaryMask(ignorePixelValueMaskGen.GetPointer()); } else { calculator->SetSecondaryMask(nullptr); } calculator->SetNBinsForHistogramStatistics(m_HistogramNBins); try { calculator->GetStatistics(); } catch (const std::exception &e) { m_LastErrorMessage = "Failure while calculating the statistics: " + std::string(e.what()); MITK_ERROR << m_LastErrorMessage; statisticCalculationSuccessful = false; } if (statisticCalculationSuccessful) { m_StatisticsContainer = calculator->GetStatistics(); auto imageRule = mitk::StatisticsToImageRelationRule::New(); imageRule->Connect(m_StatisticsContainer, m_StatisticsImage); if (nullptr != m_PlanarFigureMask) { auto maskRule = mitk::StatisticsToMaskRelationRule::New(); maskRule->Connect(m_StatisticsContainer, m_PlanarFigureMask); } if (nullptr != m_BinaryMask) { auto maskRule = mitk::StatisticsToMaskRelationRule::New(); maskRule->Connect(m_StatisticsContainer, m_BinaryMask); } m_StatisticsContainer->SetProperty(mitk::STATS_HISTOGRAM_BIN_PROPERTY_NAME.c_str(), mitk::UIntProperty::New(m_HistogramNBins)); m_StatisticsContainer->SetProperty(mitk::STATS_IGNORE_ZERO_VOXEL_PROPERTY_NAME.c_str(), mitk::BoolProperty::New(m_IgnoreZeros)); } return statisticCalculationSuccessful; } diff --git a/Modules/Segmentation/Algorithms/mitkFeatureBasedEdgeDetectionFilter.cpp b/Modules/Segmentation/Algorithms/mitkFeatureBasedEdgeDetectionFilter.cpp index ae2d461fcd..969e1c7e84 100644 --- a/Modules/Segmentation/Algorithms/mitkFeatureBasedEdgeDetectionFilter.cpp +++ b/Modules/Segmentation/Algorithms/mitkFeatureBasedEdgeDetectionFilter.cpp @@ -1,195 +1,196 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "mitkFeatureBasedEdgeDetectionFilter.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include mitk::FeatureBasedEdgeDetectionFilter::FeatureBasedEdgeDetectionFilter() { this->SetNumberOfRequiredInputs(1); this->SetNumberOfIndexedOutputs(1); } mitk::FeatureBasedEdgeDetectionFilter::~FeatureBasedEdgeDetectionFilter() { } void mitk::FeatureBasedEdgeDetectionFilter::GenerateData() { mitk::Image::ConstPointer image = ImageToUnstructuredGridFilter::GetInput(); if (m_SegmentationMask.IsNull()) { MITK_WARN << "Please set a segmentation mask first" << std::endl; return; } // First create a threshold segmentation of the image. The threshold is determined // by the mean +/- stddev of the pixel values that are covered by the segmentation mask // Compute mean and stdDev based on the current segmentation mitk::ImageStatisticsCalculator::Pointer statCalc = mitk::ImageStatisticsCalculator::New(); statCalc->SetInputImage(image); mitk::ImageMaskGenerator::Pointer imgMask = mitk::ImageMaskGenerator::New(); + imgMask->SetInputImage(image); imgMask->SetImageMask(m_SegmentationMask); auto stats = statCalc->GetStatistics()->GetStatisticsForTimeStep(0); double mean = stats.GetValueConverted(mitk::ImageStatisticsConstants::MEAN()); double stdDev = stats.GetValueConverted(mitk::ImageStatisticsConstants::STANDARDDEVIATION()); double upperThreshold = mean + stdDev; double lowerThreshold = mean - stdDev; // Perform thresholding mitk::Image::Pointer thresholdImage = mitk::Image::New(); AccessByItk_3(image.GetPointer(), ITKThresholding, lowerThreshold, upperThreshold, thresholdImage) mitk::ProgressBar::GetInstance() ->Progress(2); // Postprocess threshold segmentation // First a closing will be executed mitk::Image::Pointer closedImage = mitk::Image::New(); AccessByItk_1(thresholdImage, ThreadedClosing, closedImage); // Then we will holes that might exist mitk::MorphologicalOperations::FillHoles(closedImage); mitk::ProgressBar::GetInstance()->Progress(); // Extract the binary edges of the resulting segmentation mitk::Image::Pointer edgeImage = mitk::Image::New(); AccessByItk_1(closedImage, ContourSearch, edgeImage); // Convert the edge image into an unstructured grid mitk::ImageToUnstructuredGridFilter::Pointer i2UFilter = mitk::ImageToUnstructuredGridFilter::New(); i2UFilter->SetInput(edgeImage); i2UFilter->SetThreshold(1.0); i2UFilter->Update(); m_PointGrid = this->GetOutput(); if (m_PointGrid.IsNull()) m_PointGrid = mitk::UnstructuredGrid::New(); m_PointGrid->SetVtkUnstructuredGrid(i2UFilter->GetOutput()->GetVtkUnstructuredGrid()); mitk::ProgressBar::GetInstance()->Progress(); } template void mitk::FeatureBasedEdgeDetectionFilter::ThreadedClosing(itk::Image *originalImage, mitk::Image::Pointer &result) { typedef itk::BinaryBallStructuringElement myKernelType; myKernelType ball; ball.SetRadius(1); ball.CreateStructuringElement(); typedef typename itk::Image ImageType; typename itk::DilateObjectMorphologyImageFilter::Pointer dilationFilter = itk::DilateObjectMorphologyImageFilter::New(); dilationFilter->SetInput(originalImage); dilationFilter->SetKernel(ball); dilationFilter->Update(); typename itk::Image::Pointer dilatedImage = dilationFilter->GetOutput(); typename itk::ErodeObjectMorphologyImageFilter::Pointer erodeFilter = itk::ErodeObjectMorphologyImageFilter::New(); erodeFilter->SetInput(dilatedImage); erodeFilter->SetKernel(ball); erodeFilter->Update(); mitk::GrabItkImageMemory(erodeFilter->GetOutput(), result); } template void mitk::FeatureBasedEdgeDetectionFilter::ContourSearch(itk::Image *originalImage, mitk::Image::Pointer &result) { typedef itk::Image ImageType; typedef itk::BinaryContourImageFilter binaryContourImageFilterType; typename binaryContourImageFilterType::Pointer binaryContourFilter = binaryContourImageFilterType::New(); binaryContourFilter->SetInput(originalImage); binaryContourFilter->SetForegroundValue(1); binaryContourFilter->SetBackgroundValue(0); binaryContourFilter->Update(); typename itk::Image::Pointer itkImage = itk::Image::New(); itkImage->Graft(binaryContourFilter->GetOutput()); mitk::GrabItkImageMemory(itkImage, result); } template void mitk::FeatureBasedEdgeDetectionFilter::ITKThresholding(const itk::Image *originalImage, double lower, double upper, mitk::Image::Pointer &result) { typedef itk::Image ImageType; typedef itk::Image SegmentationType; typedef itk::BinaryThresholdImageFilter ThresholdFilterType; if (typeid(TPixel) != typeid(float) && typeid(TPixel) != typeid(double)) { // round the thresholds if we have nor a float or double image lower = std::floor(lower + 0.5); upper = std::floor(upper - 0.5); } if (lower >= upper) { upper = lower; } typename ThresholdFilterType::Pointer filter = ThresholdFilterType::New(); filter->SetInput(originalImage); filter->SetLowerThreshold(lower); filter->SetUpperThreshold(upper); filter->SetInsideValue(1); filter->SetOutsideValue(0); filter->Update(); mitk::GrabItkImageMemory(filter->GetOutput(), result); } void mitk::FeatureBasedEdgeDetectionFilter::SetSegmentationMask(mitk::Image::Pointer segmentation) { this->m_SegmentationMask = segmentation; } void mitk::FeatureBasedEdgeDetectionFilter::GenerateOutputInformation() { Superclass::GenerateOutputInformation(); }