diff --git a/Modules/Classification/CLMiniApps/CLGlobalImageFeatures.cpp b/Modules/Classification/CLMiniApps/CLGlobalImageFeatures.cpp index 3da1fcf4d7..46012dc902 100644 --- a/Modules/Classification/CLMiniApps/CLGlobalImageFeatures.cpp +++ b/Modules/Classification/CLMiniApps/CLGlobalImageFeatures.cpp @@ -1,660 +1,667 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef mitkCLPolyToNrrd_cpp #define mitkCLPolyToNrrd_cpp #include "time.h" #include #include #include #include "mitkCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "itkNearestNeighborInterpolateImageFunction.h" #include "itkResampleImageFilter.h" #include #include #include "QmitkRegisterClasses.h" #include "QmitkRenderWindow.h" #include "vtkRenderLargeImage.h" #include "vtkPNGWriter.h" typedef itk::Image< double, 3 > FloatImageType; typedef itk::Image< unsigned char, 3 > MaskImageType; template void ResampleImage(itk::Image* itkImage, float resolution, mitk::Image::Pointer& newImage) { typedef itk::Image ImageType; typedef itk::ResampleImageFilter ResampleFilterType; typename ResampleFilterType::Pointer resampler = ResampleFilterType::New(); auto spacing = itkImage->GetSpacing(); auto size = itkImage->GetLargestPossibleRegion().GetSize(); for (int i = 0; i < VImageDimension; ++i) { size[i] = size[i] / (1.0*resolution)*(1.0*spacing[i])+1.0; } spacing.Fill(resolution); resampler->SetInput(itkImage); resampler->SetSize(size); resampler->SetOutputSpacing(spacing); resampler->SetOutputOrigin(itkImage->GetOrigin()); resampler->SetOutputDirection(itkImage->GetDirection()); resampler->Update(); newImage->InitializeByItk(resampler->GetOutput()); mitk::GrabItkImageMemory(resampler->GetOutput(), newImage); } template static void CreateNoNaNMask(itk::Image* itkValue, mitk::Image::Pointer mask, mitk::Image::Pointer& newMask) { typedef itk::Image< TPixel, VImageDimension> LFloatImageType; typedef itk::Image< unsigned char, VImageDimension> LMaskImageType; typename LMaskImageType::Pointer itkMask = LMaskImageType::New(); mitk::CastToItkImage(mask, itkMask); typedef itk::ImageDuplicator< LMaskImageType > DuplicatorType; typename DuplicatorType::Pointer duplicator = DuplicatorType::New(); duplicator->SetInputImage(itkMask); duplicator->Update(); auto tmpMask = duplicator->GetOutput(); itk::ImageRegionIterator mask1Iter(itkMask, itkMask->GetLargestPossibleRegion()); itk::ImageRegionIterator mask2Iter(tmpMask, tmpMask->GetLargestPossibleRegion()); itk::ImageRegionIterator imageIter(itkValue, itkValue->GetLargestPossibleRegion()); while (!mask1Iter.IsAtEnd()) { mask2Iter.Set(0); if (mask1Iter.Value() > 0) { // Is not NaN if (imageIter.Value() == imageIter.Value()) { mask2Iter.Set(1); } } ++mask1Iter; ++mask2Iter; ++imageIter; } newMask->InitializeByItk(tmpMask); mitk::GrabItkImageMemory(tmpMask, newMask); } template static void ResampleMask(itk::Image* itkMoving, mitk::Image::Pointer ref, mitk::Image::Pointer& newMask) { typedef itk::Image< TPixel, VImageDimension> LMaskImageType; typedef itk::NearestNeighborInterpolateImageFunction< LMaskImageType> NearestNeighborInterpolateImageFunctionType; typedef itk::ResampleImageFilter ResampleFilterType; typename NearestNeighborInterpolateImageFunctionType::Pointer nn_interpolator = NearestNeighborInterpolateImageFunctionType::New(); typename LMaskImageType::Pointer itkRef = LMaskImageType::New(); mitk::CastToItkImage(ref, itkRef); typename ResampleFilterType::Pointer resampler = ResampleFilterType::New(); resampler->SetInput(itkMoving); resampler->SetReferenceImage(itkRef); resampler->UseReferenceImageOn(); resampler->SetInterpolator(nn_interpolator); resampler->Update(); newMask->InitializeByItk(resampler->GetOutput()); mitk::GrabItkImageMemory(resampler->GetOutput(), newMask); } static void ExtractSlicesFromImages(mitk::Image::Pointer image, mitk::Image::Pointer mask, mitk::Image::Pointer maskNoNaN, int direction, std::vector &imageVector, std::vector &maskVector, std::vector &maskNoNaNVector) { typedef itk::Image< double, 2 > FloatImage2DType; typedef itk::Image< unsigned char, 2 > MaskImage2DType; FloatImageType::Pointer itkFloat = FloatImageType::New(); MaskImageType::Pointer itkMask = MaskImageType::New(); MaskImageType::Pointer itkMaskNoNaN = MaskImageType::New(); mitk::CastToItkImage(mask, itkMask); mitk::CastToItkImage(maskNoNaN, itkMaskNoNaN); mitk::CastToItkImage(image, itkFloat); int idxA, idxB, idxC; switch (direction) { case 0: idxA = 1; idxB = 2; idxC = 0; break; case 1: idxA = 0; idxB = 2; idxC = 1; break; case 2: idxA = 0; idxB = 1; idxC = 2; break; default: idxA = 1; idxB = 2; idxC = 0; break; } auto imageSize = image->GetLargestPossibleRegion().GetSize(); FloatImageType::IndexType index3D; FloatImage2DType::IndexType index2D; FloatImage2DType::SpacingType spacing2D; spacing2D[0] = itkFloat->GetSpacing()[idxA]; spacing2D[1] = itkFloat->GetSpacing()[idxB]; for (int i = 0; i < imageSize[idxC]; ++i) { FloatImage2DType::RegionType region; FloatImage2DType::IndexType start; FloatImage2DType::SizeType size; start[0] = 0; start[1] = 0; size[0] = imageSize[idxA]; size[1] = imageSize[idxB]; region.SetIndex(start); region.SetSize(size); FloatImage2DType::Pointer image2D = FloatImage2DType::New(); image2D->SetRegions(region); image2D->Allocate(); MaskImage2DType::Pointer mask2D = MaskImage2DType::New(); mask2D->SetRegions(region); mask2D->Allocate(); MaskImage2DType::Pointer masnNoNaN2D = MaskImage2DType::New(); masnNoNaN2D->SetRegions(region); masnNoNaN2D->Allocate(); + unsigned long voxelsInMask = 0; + for (int a = 0; a < imageSize[idxA]; ++a) { for (int b = 0; b < imageSize[idxB]; ++b) { index3D[idxA] = a; index3D[idxB] = b; index3D[idxC] = i; index2D[0] = a; index2D[1] = b; image2D->SetPixel(index2D, itkFloat->GetPixel(index3D)); mask2D->SetPixel(index2D, itkMask->GetPixel(index3D)); masnNoNaN2D->SetPixel(index2D, itkMaskNoNaN->GetPixel(index3D)); + voxelsInMask += (itkMask->GetPixel(index3D) > 0) ? 1 : 0; } } image2D->SetSpacing(spacing2D); mask2D->SetSpacing(spacing2D); masnNoNaN2D->SetSpacing(spacing2D); mitk::Image::Pointer tmpFloatImage = mitk::Image::New(); tmpFloatImage->InitializeByItk(image2D.GetPointer()); mitk::GrabItkImageMemory(image2D, tmpFloatImage); mitk::Image::Pointer tmpMaskImage = mitk::Image::New(); tmpMaskImage->InitializeByItk(mask2D.GetPointer()); mitk::GrabItkImageMemory(mask2D, tmpMaskImage); mitk::Image::Pointer tmpMaskNoNaNImage = mitk::Image::New(); tmpMaskNoNaNImage->InitializeByItk(masnNoNaN2D.GetPointer()); mitk::GrabItkImageMemory(masnNoNaN2D, tmpMaskNoNaNImage); - imageVector.push_back(tmpFloatImage); - maskVector.push_back(tmpMaskImage); - maskNoNaNVector.push_back(tmpMaskNoNaNImage); + if (voxelsInMask > 0) + { + imageVector.push_back(tmpFloatImage); + maskVector.push_back(tmpMaskImage); + maskNoNaNVector.push_back(tmpMaskNoNaNImage); + } } } static void SaveSliceOrImageAsPNG(mitk::Image::Pointer image, mitk::Image::Pointer mask, std::string path, int index) { // Create a Standalone Datastorage for the single purpose of saving screenshots.. mitk::StandaloneDataStorage::Pointer ds = mitk::StandaloneDataStorage::New(); QmitkRenderWindow renderWindow; renderWindow.GetRenderer()->SetDataStorage(ds); auto nodeI = mitk::DataNode::New(); nodeI->SetData(image); auto nodeM = mitk::DataNode::New(); nodeM->SetData(mask); ds->Add(nodeI); ds->Add(nodeM); mitk::TimeGeometry::Pointer geo = ds->ComputeBoundingGeometry3D(ds->GetAll()); mitk::RenderingManager::GetInstance()->InitializeViews(geo); mitk::SliceNavigationController::Pointer sliceNaviController = renderWindow.GetSliceNavigationController(); unsigned int numberOfSteps = 1; if (sliceNaviController) { numberOfSteps = sliceNaviController->GetSlice()->GetSteps(); sliceNaviController->GetSlice()->SetPos(0); } renderWindow.show(); renderWindow.resize(256, 256); for (unsigned int currentStep = 0; currentStep < numberOfSteps; ++currentStep) { if (sliceNaviController) { sliceNaviController->GetSlice()->SetPos(currentStep); } renderWindow.GetRenderer()->PrepareRender(); vtkRenderWindow* renderWindow2 = renderWindow.GetVtkRenderWindow(); mitk::BaseRenderer* baserenderer = mitk::BaseRenderer::GetInstance(renderWindow2); auto vtkRender = baserenderer->GetVtkRenderer(); vtkRender->GetRenderWindow()->WaitForCompletion(); vtkRenderLargeImage* magnifier = vtkRenderLargeImage::New(); magnifier->SetInput(vtkRender); magnifier->SetMagnification(3.0); std::stringstream ss; ss << path << "_Idx-" << index << "_Step-"<> tmpImageName; auto fileWriter = vtkPNGWriter::New(); fileWriter->SetInputConnection(magnifier->GetOutputPort()); fileWriter->SetFileName(tmpImageName.c_str()); fileWriter->Write(); fileWriter->Delete(); } } int main(int argc, char* argv[]) { mitk::GIFFirstOrderStatistics::Pointer firstOrderCalculator = mitk::GIFFirstOrderStatistics::New(); mitk::GIFVolumetricStatistics::Pointer volCalculator = mitk::GIFVolumetricStatistics::New(); mitk::GIFCooccurenceMatrix::Pointer coocCalculator = mitk::GIFCooccurenceMatrix::New(); mitk::GIFCooccurenceMatrix2::Pointer cooc2Calculator = mitk::GIFCooccurenceMatrix2::New(); mitk::GIFNeighbouringGreyLevelDependenceFeature::Pointer ngldCalculator = mitk::GIFNeighbouringGreyLevelDependenceFeature::New(); mitk::GIFGrayLevelRunLength::Pointer rlCalculator = mitk::GIFGrayLevelRunLength::New(); mitk::GIFGreyLevelSizeZone::Pointer glszCalculator = mitk::GIFGreyLevelSizeZone::New(); std::vector features; features.push_back(firstOrderCalculator.GetPointer()); features.push_back(volCalculator.GetPointer()); features.push_back(coocCalculator.GetPointer()); features.push_back(cooc2Calculator.GetPointer()); features.push_back(ngldCalculator.GetPointer()); features.push_back(rlCalculator.GetPointer()); features.push_back(glszCalculator.GetPointer()); mitkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); mitk::cl::GlobalImageFeaturesParameter param; param.AddParameter(parser); parser.addArgument("--","-", mitkCommandLineParser::String, "---", "---", us::Any(),true); for (auto cFeature : features) { cFeature->AddArguments(parser); } parser.addArgument("--", "-", mitkCommandLineParser::String, "---", "---", us::Any(), true); parser.addArgument("description","d",mitkCommandLineParser::String,"Text","Description that is added to the output",us::Any()); parser.addArgument("direction", "dir", mitkCommandLineParser::String, "Int", "Allows to specify the direction for Cooc and RL. 0: All directions, 1: Only single direction (Test purpose), 2,3,4... Without dimension 0,1,2... ", us::Any()); parser.addArgument("slice-wise", "slice", mitkCommandLineParser::String, "Int", "Allows to specify if the image is processed slice-wise (number giving direction) ", us::Any()); parser.addArgument("output-mode", "omode", mitkCommandLineParser::Int, "Int", "Defines if the results of an image / slice are written in a single row (0 , default) or column (1)."); // Miniapp Infos parser.setCategory("Classification Tools"); parser.setTitle("Global Image Feature calculator"); parser.setDescription("Calculates different global statistics for a given segmentation / image combination"); parser.setContributor("MBI"); std::map parsedArgs = parser.parseArguments(argc, argv); param.ParseParameter(parsedArgs); if (parsedArgs.size()==0) { return EXIT_FAILURE; } if ( parsedArgs.count("help") || parsedArgs.count("h")) { return EXIT_SUCCESS; } bool savePNGofSlices = true; std::string folderForPNGOfSlices = "E:\\tmp\\bonekamp\\fig\\"; std::string version = "Version: 1.19"; MITK_INFO << version; std::ofstream log; if (param.useLogfile) { log.open(param.logfilePath, std::ios::app); log << version; log << "Image: " << param.imagePath; log << "Mask: " << param.maskPath; } mitk::Image::Pointer image; mitk::Image::Pointer mask; mitk::Image::Pointer tmpImage = mitk::IOUtil::LoadImage(param.imagePath); mitk::Image::Pointer tmpMask = mitk::IOUtil::LoadImage(param.maskPath); image = tmpImage; mask = tmpMask; if ((image->GetDimension() != mask->GetDimension())) { MITK_INFO << "Dimension of image does not match. "; MITK_INFO << "Correct one image, may affect the result"; if (image->GetDimension() == 2) { mitk::Convert2Dto3DImageFilter::Pointer multiFilter2 = mitk::Convert2Dto3DImageFilter::New(); multiFilter2->SetInput(tmpImage); multiFilter2->Update(); image = multiFilter2->GetOutput(); } if (mask->GetDimension() == 2) { mitk::Convert2Dto3DImageFilter::Pointer multiFilter3 = mitk::Convert2Dto3DImageFilter::New(); multiFilter3->SetInput(tmpMask); multiFilter3->Update(); mask = multiFilter3->GetOutput(); } } int writeDirection = 0; if (parsedArgs.count("output-mode")) { writeDirection = us::any_cast(parsedArgs["output-mode"]); } if (param.resampleToFixIsotropic) { mitk::Image::Pointer newImage = mitk::Image::New(); AccessByItk_2(image, ResampleImage, param.resampleResolution, newImage); image = newImage; } if ( ! mitk::Equal(mask->GetGeometry(0)->GetOrigin(), image->GetGeometry(0)->GetOrigin())) { MITK_INFO << "Not equal Origins"; if (param.ensureSameSpace) { MITK_INFO << "Warning!"; MITK_INFO << "The origin of the input image and the mask do not match. They are"; MITK_INFO << "now corrected. Please check to make sure that the images still match"; image->GetGeometry(0)->SetOrigin(mask->GetGeometry(0)->GetOrigin()); } else { return -1; } } if (param.resampleMask) { mitk::Image::Pointer newMaskImage = mitk::Image::New(); AccessByItk_2(mask, ResampleMask, image, newMaskImage); mask = newMaskImage; } if ( ! mitk::Equal(mask->GetGeometry(0)->GetSpacing(), image->GetGeometry(0)->GetSpacing())) { MITK_INFO << "Not equal Sapcings"; if (param.ensureSameSpace) { MITK_INFO << "Warning!"; MITK_INFO << "The spacing of the mask was set to match the spacing of the input image."; MITK_INFO << "This might cause unintended spacing of the mask image"; image->GetGeometry(0)->SetSpacing(mask->GetGeometry(0)->GetSpacing()); } else { MITK_INFO << "The spacing of the mask and the input images is not equal."; MITK_INFO << "Terminating the programm. You may use the '-fi' option"; return -1; } } int direction = 0; if (parsedArgs.count("direction")) { direction = mitk::cl::splitDouble(parsedArgs["direction"].ToString(), ';')[0]; } MITK_INFO << "Start creating Mask without NaN"; mitk::Image::Pointer maskNoNaN = mitk::Image::New(); AccessByItk_2(image, CreateNoNaNMask, mask, maskNoNaN); //CreateNoNaNMask(mask, image, maskNoNaN); bool sliceWise = false; int sliceDirection = 0; int currentSlice = 0; bool imageToProcess = true; std::vector floatVector; std::vector maskVector; std::vector maskNoNaNVector; if ((parsedArgs.count("slice-wise")) && image->GetDimension() > 2) { MITK_INFO << "Enabled slice-wise"; sliceWise = true; sliceDirection = mitk::cl::splitDouble(parsedArgs["slice-wise"].ToString(), ';')[0]; MITK_INFO << sliceDirection; ExtractSlicesFromImages(image, mask, maskNoNaN, sliceDirection, floatVector, maskVector, maskNoNaNVector); MITK_INFO << "Slice"; } for (auto cFeature : features) { if (param.defineGlobalMinimumIntensity) { cFeature->SetMinimumIntensity(param.globalMinimumIntensity); cFeature->SetUseMinimumIntensity(true); } if (param.defineGlobalMaximumIntensity) { cFeature->SetMaximumIntensity(param.globalMaximumIntensity); cFeature->SetUseMaximumIntensity(true); } if (param.defineGlobalNumberOfBins) { cFeature->SetBins(param.globalNumberOfBins); + MITK_INFO << param.globalNumberOfBins; } cFeature->SetParameter(parsedArgs); cFeature->SetDirection(direction); } bool addDescription = parsedArgs.count("description"); mitk::cl::FeatureResultWritter writer(param.outputPath, writeDirection); std::string description = ""; if (addDescription) { description = parsedArgs["description"].ToString(); } mitk::Image::Pointer cImage = image; mitk::Image::Pointer cMask = mask; mitk::Image::Pointer cMaskNoNaN = maskNoNaN; if (param.useHeader) { writer.AddColumn("Patient"); writer.AddColumn("Image"); writer.AddColumn("Segmentation"); } // Create a QTApplication and a Datastorage // This is necessary in order to save screenshots of // each image / slice. QApplication qtapplication(argc, argv); QmitkRegisterClasses(); std::vector allStats; while (imageToProcess) { if (sliceWise) { cImage = floatVector[currentSlice]; cMask = maskVector[currentSlice]; cMaskNoNaN = maskNoNaNVector[currentSlice]; imageToProcess = (floatVector.size()-1 > (currentSlice)) ? true : false ; } else { imageToProcess = false; } if (param.writePNGScreenshots) { SaveSliceOrImageAsPNG(cImage, cMask, param.pngScreenshotsPath, currentSlice); } if (param.writeAnalysisImage) { mitk::IOUtil::SaveImage(cImage, param.anaylsisImagePath); } if (param.writeAnalysisMask) { mitk::IOUtil::SaveImage(cImage, param.analysisMaskPath); } mitk::AbstractGlobalImageFeature::FeatureListType stats; for (auto cFeature : features) { cFeature->CalculateFeaturesUsingParameters(cImage, cMask, cMaskNoNaN, stats); } for (std::size_t i = 0; i < stats.size(); ++i) { std::cout << stats[i].first << " - " << stats[i].second << std::endl; } writer.AddHeader(description, currentSlice, stats, param.useHeader, addDescription); if (true) { writer.AddColumn(param.imageFolder); writer.AddColumn(param.imageName); writer.AddColumn(param.maskName); } writer.AddResult(description, currentSlice, stats, param.useHeader, addDescription); allStats.push_back(stats); ++currentSlice; } if (sliceWise) { mitk::AbstractGlobalImageFeature::FeatureListType statMean, statStd; for (std::size_t i = 0; i < allStats[0].size(); ++i) { auto cElement1 = allStats[0][i]; cElement1.first = "SliceWise Mean " + cElement1.first; cElement1.second = 0.0; auto cElement2 = allStats[0][i]; cElement2.first = "SliceWise Var. " + cElement2.first; cElement2.second = 0.0; statMean.push_back(cElement1); statStd.push_back(cElement2); } for (auto cStat : allStats) { for (std::size_t i = 0; i < cStat.size(); ++i) { statMean[i].second += cStat[i].second / (1.0*allStats.size()); } } for (auto cStat : allStats) { for (std::size_t i = 0; i < cStat.size(); ++i) { statStd[i].second += (cStat[i].second - statMean[i].second)*(cStat[i].second - statMean[i].second) / (1.0*allStats.size()); } } for (std::size_t i = 0; i < statMean.size(); ++i) { std::cout << statMean[i].first << " - " << statMean[i].second << std::endl; std::cout << statStd[i].first << " - " << statStd[i].second << std::endl; } if (true) { writer.AddColumn(param.imageFolder); writer.AddColumn(param.imageName); writer.AddColumn(param.maskName + " - Mean"); } writer.AddResult(description, currentSlice, statMean, param.useHeader, addDescription); if (true) { writer.AddColumn(param.imageFolder); writer.AddColumn(param.imageName); writer.AddColumn(param.maskName + " - Var."); } writer.AddResult(description, currentSlice, statStd, param.useHeader, addDescription); } if (param.useLogfile) { log << "Finished calculation"; log.close(); } return 0; } #endif diff --git a/Modules/Classification/CLMiniApps/CLScreenshot.cpp b/Modules/Classification/CLMiniApps/CLScreenshot.cpp new file mode 100644 index 0000000000..8e51429257 --- /dev/null +++ b/Modules/Classification/CLMiniApps/CLScreenshot.cpp @@ -0,0 +1,145 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ +#ifndef mitkCLPolyToNrrd_cpp +#define mitkCLPolyToNrrd_cpp + +#include "time.h" +#include +#include + +#include +#include "mitkCommandLineParser.h" + +#include + +#include +#include +#include "QmitkRegisterClasses.h" +#include "QmitkRenderWindow.h" +#include "vtkRenderLargeImage.h" +#include "vtkPNGWriter.h" + + +static +void SaveSliceOrImageAsPNG(std::vector listOfOutputs, std::string path) +{ + // Create a Standalone Datastorage for the single purpose of saving screenshots.. + mitk::StandaloneDataStorage::Pointer ds = mitk::StandaloneDataStorage::New(); + QmitkRenderWindow renderWindow; + renderWindow.GetRenderer()->SetDataStorage(ds); + + for (auto name : listOfOutputs) + { + mitk::Image::Pointer tmpImage = mitk::IOUtil::LoadImage(name); + auto nodeI = mitk::DataNode::New(); + nodeI->SetData(tmpImage); + ds->Add(nodeI); + } + + mitk::TimeGeometry::Pointer geo = ds->ComputeBoundingGeometry3D(ds->GetAll()); + mitk::RenderingManager::GetInstance()->InitializeViews(geo); + + mitk::SliceNavigationController::Pointer sliceNaviController = renderWindow.GetSliceNavigationController(); + unsigned int numberOfSteps = 1; + if (sliceNaviController) + { + numberOfSteps = sliceNaviController->GetSlice()->GetSteps(); + sliceNaviController->GetSlice()->SetPos(0); + } + + renderWindow.show(); + renderWindow.resize(512, 512); + + for (unsigned int currentStep = 0; currentStep < numberOfSteps; ++currentStep) + { + if (sliceNaviController) + { + sliceNaviController->GetSlice()->SetPos(currentStep); + } + + renderWindow.GetRenderer()->PrepareRender(); + + vtkRenderWindow* renderWindow2 = renderWindow.GetVtkRenderWindow(); + mitk::BaseRenderer* baserenderer = mitk::BaseRenderer::GetInstance(renderWindow2); + auto vtkRender = baserenderer->GetVtkRenderer(); + vtkRender->GetRenderWindow()->WaitForCompletion(); + + vtkRenderLargeImage* magnifier = vtkRenderLargeImage::New(); + magnifier->SetInput(vtkRender); + magnifier->SetMagnification(3.0); + + std::stringstream ss; + ss << path << "screenshot_step-"<> tmpImageName; + auto fileWriter = vtkPNGWriter::New(); + fileWriter->SetInputConnection(magnifier->GetOutputPort()); + fileWriter->SetFileName(tmpImageName.c_str()); + fileWriter->Write(); + fileWriter->Delete(); + } +} + +int main(int argc, char* argv[]) +{ + mitkCommandLineParser parser; + parser.setArgumentPrefix("--", "-"); + + // Required Parameter + parser.addArgument("image", "i", mitkCommandLineParser::InputImage, "Input Image", "Path to the input image files (Separated with semicolons)", us::Any(), false); + parser.addArgument("output", "o", mitkCommandLineParser::OutputFile, "Output text file", "Path to output file. The output statistic is appended to this file.", us::Any(), false); + parser.addArgument("direction", "dir", mitkCommandLineParser::String, "Int", "Allows to specify the direction for Cooc and RL. 0: All directions, 1: Only single direction (Test purpose), 2,3,4... Without dimension 0,1,2... ", us::Any()); + + // Miniapp Infos + parser.setCategory("Classification Tools"); + parser.setTitle("Screenshot of a single image"); + parser.setDescription(""); + parser.setContributor("MBI"); + + std::map parsedArgs = parser.parseArguments(argc, argv); + + if (parsedArgs.size()==0) + { + return EXIT_FAILURE; + } + if ( parsedArgs.count("help") || parsedArgs.count("h")) + { + return EXIT_SUCCESS; + } + + std::string version = "Version: 1.0"; + MITK_INFO << version; + + int direction = 0; + if (parsedArgs.count("direction")) + { + direction = mitk::cl::splitDouble(parsedArgs["direction"].ToString(), ';')[0]; + } + + auto listOfFiles = mitk::cl::splitString(parsedArgs["image"].ToString(), ';'); + + // Create a QTApplication and a Datastorage + // This is necessary in order to save screenshots of + // each image / slice. + QApplication qtapplication(argc, argv); + QmitkRegisterClasses(); + + SaveSliceOrImageAsPNG(listOfFiles, parsedArgs["output"].ToString()); + + return 0; +} + +#endif diff --git a/Modules/Classification/CLMiniApps/CMakeLists.txt b/Modules/Classification/CLMiniApps/CMakeLists.txt index b1010950c0..9dac692b07 100644 --- a/Modules/Classification/CLMiniApps/CMakeLists.txt +++ b/Modules/Classification/CLMiniApps/CMakeLists.txt @@ -1,139 +1,140 @@ option(BUILD_ClassificationMiniApps "Build commandline tools for classification" OFF) if(BUILD_ClassificationMiniApps OR MITK_BUILD_ALL_APPS) include_directories( ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ) # list of miniapps # if an app requires additional dependencies # they are added after a "^^" and separated by "_" set( classificationminiapps RandomForestTraining^^MitkCLVigraRandomForest NativeHeadCTSegmentation^^MitkCLVigraRandomForest ManualSegmentationEvaluation^^MitkCLVigraRandomForest + CLScreenshot^^MitkCore_MitkQtWidgetsExt_MitkCLUtilities CLGlobalImageFeatures^^MitkCore_MitkCLUtilities_MitkQtWidgetsExt CLMRNormalization^^MitkCore_MitkCLUtilities_MitkCLMRUtilities CLStaple^^MitkCore_MitkCLUtilities CLVoxelFeatures^^MitkCore_MitkCLUtilities CLDicom2Nrrd^^MitkCore CLPolyToNrrd^^MitkCore CLImageTypeConverter^^MitkCore CLResampleImageToReference^^MitkCore CLRandomSampling^^MitkCore_MitkCLUtilities CLRemoveEmptyVoxels^^MitkCore CLN4^^MitkCore CLMultiForestPrediction^^MitkDataCollection_MitkCLVigraRandomForest CLNrrdToPoly^^MitkCore CL2Dto3DImage^^MitkCore CLWeighting^^MitkCore_MitkCLImportanceWeighting_MitkCLUtilities CLOverlayRoiCenterOfMass^^MitkCore_MitkCLUtilities_MitkQtWidgetsExt # RandomForestPrediction^^MitkCLVigraRandomForest ) foreach(classificationminiapps ${classificationminiapps}) # extract mini app name and dependencies string(REPLACE "^^" "\\;" miniapp_info ${classificationminiapps}) set(miniapp_info_list ${miniapp_info}) list(GET miniapp_info_list 0 appname) list(GET miniapp_info_list 1 raw_dependencies) string(REPLACE "_" "\\;" dependencies "${raw_dependencies}") set(dependencies_list ${dependencies}) mitk_create_executable(${appname} DEPENDS MitkCore MitkCLCore MitkCommandLine ${dependencies_list} PACKAGE_DEPENDS ITK Qt5|Core Vigra CPP_FILES ${appname}.cpp ) # CPP_FILES ${appname}.cpp mitkCommandLineParser.cpp if(EXECUTABLE_IS_ENABLED) # On Linux, create a shell script to start a relocatable application if(UNIX AND NOT APPLE) install(PROGRAMS "${MITK_SOURCE_DIR}/CMake/RunInstalledApp.sh" DESTINATION "." RENAME ${EXECUTABLE_TARGET}.sh) endif() get_target_property(_is_bundle ${EXECUTABLE_TARGET} MACOSX_BUNDLE) if(APPLE) if(_is_bundle) set(_target_locations ${EXECUTABLE_TARGET}.app) set(${_target_locations}_qt_plugins_install_dir ${EXECUTABLE_TARGET}.app/Contents/MacOS) set(_bundle_dest_dir ${EXECUTABLE_TARGET}.app/Contents/MacOS) set(_qt_plugins_for_current_bundle ${EXECUTABLE_TARGET}.app/Contents/MacOS) set(_qt_conf_install_dirs ${EXECUTABLE_TARGET}.app/Contents/Resources) install(TARGETS ${EXECUTABLE_TARGET} BUNDLE DESTINATION . ) else() if(NOT MACOSX_BUNDLE_NAMES) set(_qt_conf_install_dirs bin) set(_target_locations bin/${EXECUTABLE_TARGET}) set(${_target_locations}_qt_plugins_install_dir bin) install(TARGETS ${EXECUTABLE_TARGET} RUNTIME DESTINATION bin) else() foreach(bundle_name ${MACOSX_BUNDLE_NAMES}) list(APPEND _qt_conf_install_dirs ${bundle_name}.app/Contents/Resources) set(_current_target_location ${bundle_name}.app/Contents/MacOS/${EXECUTABLE_TARGET}) list(APPEND _target_locations ${_current_target_location}) set(${_current_target_location}_qt_plugins_install_dir ${bundle_name}.app/Contents/MacOS) message( " set(${_current_target_location}_qt_plugins_install_dir ${bundle_name}.app/Contents/MacOS) ") install(TARGETS ${EXECUTABLE_TARGET} RUNTIME DESTINATION ${bundle_name}.app/Contents/MacOS/) endforeach() endif() endif() else() set(_target_locations bin/${EXECUTABLE_TARGET}${CMAKE_EXECUTABLE_SUFFIX}) set(${_target_locations}_qt_plugins_install_dir bin) set(_qt_conf_install_dirs bin) install(TARGETS ${EXECUTABLE_TARGET} RUNTIME DESTINATION bin) endif() endif() endforeach() # This mini app does not depend on mitkDiffusionImaging at all #mitk_create_executable(CLGlobalImageFeatures # DEPENDS MitkCore MitkCLUtilities # CPP_FILES CLGlobalImageFeatures.cpp mitkCommandLineParser.cpp #) mitk_create_executable(CLSimpleVoxelClassification DEPENDS MitkCore MitkCLCore MitkDataCollection MitkCLVigraRandomForest MitkCommandLine CPP_FILES CLSimpleVoxelClassification.cpp ) # This mini app does not depend on mitkDiffusionImaging at all mitk_create_executable(CLVoxelClassification DEPENDS MitkCore MitkCLCore MitkDataCollection MitkCLImportanceWeighting MitkCLVigraRandomForest CPP_FILES CLVoxelClassification.cpp ) #mitk_create_executable(CLBrainMask # DEPENDS MitkCore MitkCLUtilities # CPP_FILES CLBrainMask.cpp mitkCommandLineParser.cpp #) #mitk_create_executable(CLImageConverter # DEPENDS MitkCore # CPP_FILES CLImageConverter.cpp mitkCommandLineParser.cpp #) #mitk_create_executable(CLSurWeighting # DEPENDS MitkCore MitkCLUtilities MitkDataCollection #MitkCLImportanceWeighting # CPP_FILES CLSurWeighting.cpp mitkCommandLineParser.cpp #) #mitk_create_executable(CLImageCropper # DEPENDS MitkCore MitkCLUtilities MitkAlgorithmsExt # CPP_FILES CLImageCropper.cpp mitkCommandLineParser.cpp #) # On Linux, create a shell script to start a relocatable application if(UNIX AND NOT APPLE) install(PROGRAMS "${MITK_SOURCE_DIR}/CMake/RunInstalledApp.sh" DESTINATION "." RENAME ${EXECUTABLE_TARGET}.sh) endif() if(EXECUTABLE_IS_ENABLED) MITK_INSTALL_TARGETS(EXECUTABLES ${EXECUTABLE_TARGET}) endif() endif() diff --git a/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToRunLengthFeaturesFilter.hxx b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToRunLengthFeaturesFilter.hxx index 3b8d3d5ce6..c38efff322 100644 --- a/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToRunLengthFeaturesFilter.hxx +++ b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToRunLengthFeaturesFilter.hxx @@ -1,663 +1,663 @@ /*=================================================================== 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. ===================================================================*/ /*========================================================================= * * Copyright Insight Software Consortium * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0.txt * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. * *=========================================================================*/ #ifndef __itkEnhancedHistogramToRunLengthFeaturesFilter_hxx #define __itkEnhancedHistogramToRunLengthFeaturesFilter_hxx #include "itkEnhancedHistogramToRunLengthFeaturesFilter.h" #include "itkNumericTraits.h" #include "vnl/vnl_math.h" namespace itk { namespace Statistics { //constructor template EnhancedHistogramToRunLengthFeaturesFilter ::EnhancedHistogramToRunLengthFeaturesFilter() : m_NumberOfVoxels(1) { this->ProcessObject::SetNumberOfRequiredInputs( 1 ); // allocate the data objects for the outputs which are // just decorators real types for( unsigned int i = 0; i < 17; i++ ) { this->ProcessObject::SetNthOutput( i, this->MakeOutput( i ) ); } } template void EnhancedHistogramToRunLengthFeaturesFilter< THistogram> ::SetInput( const HistogramType *histogram ) { this->ProcessObject::SetNthInput( 0, const_cast( histogram ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::HistogramType * EnhancedHistogramToRunLengthFeaturesFilter< THistogram> ::GetInput() const { if ( this->GetNumberOfInputs() < 1 ) { return ITK_NULLPTR; } return itkDynamicCastInDebugMode(this->ProcessObject::GetInput( 0 ) ); } template typename EnhancedHistogramToRunLengthFeaturesFilter::DataObjectPointer EnhancedHistogramToRunLengthFeaturesFilter ::MakeOutput( DataObjectPointerArraySizeType itkNotUsed( idx ) ) { return MeasurementObjectType::New().GetPointer(); } template void EnhancedHistogramToRunLengthFeaturesFilter< THistogram>:: GenerateData( void ) { const HistogramType * inputHistogram = this->GetInput(); this->m_TotalNumberOfRuns = static_cast ( inputHistogram->GetTotalFrequency() ); MeasurementType shortRunEmphasis = NumericTraits::ZeroValue(); MeasurementType longRunEmphasis = NumericTraits::ZeroValue(); MeasurementType greyLevelNonuniformity = NumericTraits::ZeroValue(); MeasurementType runLengthNonuniformity = NumericTraits::ZeroValue(); MeasurementType lowGreyLevelRunEmphasis = NumericTraits::ZeroValue(); MeasurementType highGreyLevelRunEmphasis = NumericTraits::ZeroValue(); MeasurementType shortRunLowGreyLevelEmphasis = NumericTraits::ZeroValue(); MeasurementType shortRunHighGreyLevelEmphasis = NumericTraits::ZeroValue(); MeasurementType longRunLowGreyLevelEmphasis = NumericTraits::ZeroValue(); MeasurementType longRunHighGreyLevelEmphasis = NumericTraits::ZeroValue(); MeasurementType runPercentage = NumericTraits::ZeroValue(); MeasurementType numberOfRuns = NumericTraits::ZeroValue(); //Added 15.07.2016 MeasurementType greyLevelVariance = NumericTraits::ZeroValue(); MeasurementType runLengthVariance = NumericTraits::ZeroValue(); MeasurementType runEntropy = NumericTraits::ZeroValue(); //Added 09.09.2016 MeasurementType greyLevelNonuniformityNormalized = NumericTraits::ZeroValue(); MeasurementType runLengthNonuniformityNormalized = NumericTraits::ZeroValue(); vnl_vector greyLevelNonuniformityVector( inputHistogram->GetSize()[0], 0.0 ); vnl_vector runLengthNonuniformityVector( inputHistogram->GetSize()[1], 0.0 ); typedef typename HistogramType::ConstIterator HistogramIterator; double mu_i = 0.0; double mu_j = 0.0; //Calculate the means. for ( HistogramIterator hit = inputHistogram->Begin(); hit != inputHistogram->End(); ++hit ) { IndexType index = hit.GetIndex(); MeasurementType frequency = hit.GetFrequency(); if ( frequency == 0 ) { continue; } // MITK_INFO << index[0] + 1 << "|" << index[1] + 1 << " " << frequency; MeasurementVectorType measurement = hit.GetMeasurementVector(); double i = index[0] + 1; double j = index[1] + 1; double p_ij = frequency / (1.0*m_TotalNumberOfRuns); mu_i += i * p_ij; mu_j += j * p_ij; } // MITK_INFO << "Mu_i " << mu_i << " Mu_j " << mu_j; //Calculate the other features. const double log2 = std::log(2.0); for ( HistogramIterator hit = inputHistogram->Begin(); hit != inputHistogram->End(); ++hit ) { MeasurementType frequency = hit.GetFrequency(); if ( frequency == 0 ) { continue; } MeasurementVectorType measurement = hit.GetMeasurementVector(); IndexType index = hit.GetIndex(); // inputHistogram->GetIndex( hit.GetInstanceIdentifier() ); - double i2 = static_cast( ( index[0] + 1 ) * ( index[0] + 1 ) ); - double j2 = static_cast( ( index[1] + 1 ) * ( index[1] + 1 ) ); - double i = index[0] + 1; double j = index[1] + 1; + double i2 = i*i; + double j2 = j*j; + MITK_INFO << i << " " << j; double p_ij = frequency / m_TotalNumberOfRuns; greyLevelVariance += ((i - mu_i) * (i - mu_i) * p_ij); runLengthVariance += ((j - mu_j) * (j - mu_j) * p_ij); runEntropy -= ( p_ij > 0.0001 ) ? p_ij *std::log(p_ij) / log2 : 0; // Traditional measures shortRunEmphasis += ( frequency / j2 ); longRunEmphasis += ( frequency * j2 ); greyLevelNonuniformityVector[index[0]] += frequency; runLengthNonuniformityVector[index[1]] += frequency; // measures from Chu et al. lowGreyLevelRunEmphasis += (i2 > 0.0001) ? ( frequency / i2 ) : 0; highGreyLevelRunEmphasis += ( frequency * i2 ); // measures from Dasarathy and Holder shortRunLowGreyLevelEmphasis += ((i2 * j2) > 0.0001) ? ( frequency / ( i2 * j2 ) ) : 0; shortRunHighGreyLevelEmphasis += (j2 > 0.0001) ? ( frequency * i2 / j2 ) : 0; longRunLowGreyLevelEmphasis += (i2 > 0.0001) ? ( frequency * j2 / i2 ) : 0; longRunHighGreyLevelEmphasis += ( frequency * i2 * j2 ); } greyLevelNonuniformity = greyLevelNonuniformityVector.squared_magnitude(); runLengthNonuniformity = runLengthNonuniformityVector.squared_magnitude(); // Normalize all measures by the total number of runs if (this->m_TotalNumberOfRuns > 0) { shortRunEmphasis /= static_cast( this->m_TotalNumberOfRuns ); longRunEmphasis /= static_cast( this->m_TotalNumberOfRuns ); greyLevelNonuniformity /= static_cast( this->m_TotalNumberOfRuns ); runLengthNonuniformity /= static_cast( this->m_TotalNumberOfRuns ); lowGreyLevelRunEmphasis /= static_cast( this->m_TotalNumberOfRuns ); highGreyLevelRunEmphasis /= static_cast( this->m_TotalNumberOfRuns ); shortRunLowGreyLevelEmphasis /= static_cast( this->m_TotalNumberOfRuns ); shortRunHighGreyLevelEmphasis /= static_cast( this->m_TotalNumberOfRuns ); longRunLowGreyLevelEmphasis /= static_cast( this->m_TotalNumberOfRuns ); longRunHighGreyLevelEmphasis /= static_cast( this->m_TotalNumberOfRuns ); runPercentage = static_cast( this->m_TotalNumberOfRuns ) / static_cast( this->m_NumberOfVoxels ); numberOfRuns = static_cast( this->m_TotalNumberOfRuns ) ; greyLevelNonuniformityNormalized = greyLevelNonuniformity / static_cast(this->m_TotalNumberOfRuns); runLengthNonuniformityNormalized = runLengthNonuniformity / static_cast(this->m_TotalNumberOfRuns); } else { shortRunEmphasis = 0; longRunEmphasis = 0; greyLevelNonuniformity = 0; runLengthNonuniformity= 0; lowGreyLevelRunEmphasis = 0; highGreyLevelRunEmphasis = 0; shortRunLowGreyLevelEmphasis = 0; shortRunHighGreyLevelEmphasis= 0; longRunLowGreyLevelEmphasis = 0; longRunHighGreyLevelEmphasis = 0; runPercentage = 0; greyLevelNonuniformityNormalized = 0; runLengthNonuniformityNormalized = 0; numberOfRuns = static_cast( this->m_TotalNumberOfRuns ) ; } MeasurementObjectType* shortRunEmphasisOutputObject = static_cast( this->ProcessObject::GetOutput( 0 ) ); shortRunEmphasisOutputObject->Set( shortRunEmphasis ); MeasurementObjectType* longRunEmphasisOutputObject = static_cast( this->ProcessObject::GetOutput( 1 ) ); longRunEmphasisOutputObject->Set( longRunEmphasis ); MeasurementObjectType* greyLevelNonuniformityOutputObject = static_cast( this->ProcessObject::GetOutput( 2 ) ); greyLevelNonuniformityOutputObject->Set( greyLevelNonuniformity ); MeasurementObjectType* greyLevelNonuniformityNormalizedOutputObject = static_cast( this->ProcessObject::GetOutput( 15 ) ); greyLevelNonuniformityNormalizedOutputObject->Set( greyLevelNonuniformityNormalized ); MeasurementObjectType* runLengthNonuniformityNormalizedOutputObject = static_cast( this->ProcessObject::GetOutput( 16 ) ); runLengthNonuniformityNormalizedOutputObject->Set( runLengthNonuniformityNormalized ); MeasurementObjectType* runLengthNonuniformityOutputObject = static_cast( this->ProcessObject::GetOutput( 3 ) ); runLengthNonuniformityOutputObject->Set( runLengthNonuniformity ); MeasurementObjectType* lowGreyLevelRunEmphasisOutputObject = static_cast( this->ProcessObject::GetOutput( 4 ) ); lowGreyLevelRunEmphasisOutputObject->Set( lowGreyLevelRunEmphasis ); MeasurementObjectType* highGreyLevelRunEmphasisOutputObject = static_cast( this->ProcessObject::GetOutput( 5 ) ); highGreyLevelRunEmphasisOutputObject->Set( highGreyLevelRunEmphasis ); MeasurementObjectType* shortRunLowGreyLevelEmphasisOutputObject = static_cast( this->ProcessObject::GetOutput( 6 ) ); shortRunLowGreyLevelEmphasisOutputObject->Set( shortRunLowGreyLevelEmphasis ); MeasurementObjectType* shortRunHighGreyLevelEmphasisOutputObject = static_cast( this->ProcessObject::GetOutput( 7 ) ); shortRunHighGreyLevelEmphasisOutputObject->Set( shortRunHighGreyLevelEmphasis ); MeasurementObjectType* longRunLowGreyLevelEmphasisOutputObject = static_cast( this->ProcessObject::GetOutput( 8 ) ); longRunLowGreyLevelEmphasisOutputObject->Set( longRunLowGreyLevelEmphasis ); MeasurementObjectType* longRunHighGreyLevelEmphasisOutputObject = static_cast( this->ProcessObject::GetOutput( 9 ) ); longRunHighGreyLevelEmphasisOutputObject->Set( longRunHighGreyLevelEmphasis ); MeasurementObjectType* runPercentagesOutputObject = static_cast( this->ProcessObject::GetOutput( 10 ) ); runPercentagesOutputObject->Set( runPercentage ); MeasurementObjectType* numberOfRunsOutputObject = static_cast( this->ProcessObject::GetOutput( 11 ) ); numberOfRunsOutputObject->Set( numberOfRuns ); MeasurementObjectType* greyLevelVarianceOutputObject = static_cast( this->ProcessObject::GetOutput( 12 ) ); greyLevelVarianceOutputObject->Set( greyLevelVariance ); MeasurementObjectType* runLengthVarianceOutputObject = static_cast( this->ProcessObject::GetOutput( 13 ) ); runLengthVarianceOutputObject->Set( runLengthVariance ); MeasurementObjectType* runEntropyOutputObject = static_cast( this->ProcessObject::GetOutput( 14 ) ); runEntropyOutputObject->Set( runEntropy ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetShortRunEmphasisOutput() const { return itkDynamicCastInDebugMode(this->ProcessObject::GetOutput( 0 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetLongRunEmphasisOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 1 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetGreyLevelNonuniformityOutput() const { return itkDynamicCastInDebugMode(this->ProcessObject::GetOutput( 2 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetRunLengthNonuniformityOutput() const { return itkDynamicCastInDebugMode(this->ProcessObject::GetOutput( 3 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetLowGreyLevelRunEmphasisOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 4 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetHighGreyLevelRunEmphasisOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 5 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetShortRunLowGreyLevelEmphasisOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 6 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetShortRunHighGreyLevelEmphasisOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 7 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetLongRunLowGreyLevelEmphasisOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 8 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetLongRunHighGreyLevelEmphasisOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 9 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetRunPercentageOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 10 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetNumberOfRunsOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 11 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetGreyLevelVarianceOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 12 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetRunLengthVarianceOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 13 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetRunEntropyOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 14 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetGreyLevelNonuniformityNormalizedOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 15 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetRunLengthNonuniformityNormalizedOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 16 ) ); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetShortRunEmphasis() const { return this->GetShortRunEmphasisOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetLongRunEmphasis() const { return this->GetLongRunEmphasisOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetGreyLevelNonuniformity() const { return this->GetGreyLevelNonuniformityOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetRunLengthNonuniformity() const { return this->GetRunLengthNonuniformityOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetLowGreyLevelRunEmphasis() const { return this->GetLowGreyLevelRunEmphasisOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetHighGreyLevelRunEmphasis() const { return this->GetHighGreyLevelRunEmphasisOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetShortRunLowGreyLevelEmphasis() const { return this->GetShortRunLowGreyLevelEmphasisOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetShortRunHighGreyLevelEmphasis() const { return this->GetShortRunHighGreyLevelEmphasisOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetLongRunLowGreyLevelEmphasis() const { return this->GetLongRunLowGreyLevelEmphasisOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetLongRunHighGreyLevelEmphasis() const { return this->GetLongRunHighGreyLevelEmphasisOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetRunPercentage() const { return this->GetRunPercentageOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetNumberOfRuns() const { return this->GetNumberOfRunsOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetGreyLevelVariance() const { return this->GetGreyLevelVarianceOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetRunLengthVariance() const { return this->GetRunLengthVarianceOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetRunEntropy() const { return this->GetRunEntropyOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetGreyLevelNonuniformityNormalized() const { return this->GetGreyLevelNonuniformityNormalizedOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetRunLengthNonuniformityNormalized() const { return this->GetRunLengthNonuniformityNormalizedOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter< THistogram>::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetFeature( RunLengthFeatureName feature ) { switch( feature ) { case ShortRunEmphasis: return this->GetShortRunEmphasis(); case LongRunEmphasis: return this->GetLongRunEmphasis(); case GreyLevelNonuniformity: return this->GetGreyLevelNonuniformity(); case GreyLevelNonuniformityNormalized: return this->GetGreyLevelNonuniformityNormalized(); case RunLengthNonuniformity: return this->GetRunLengthNonuniformity(); case RunLengthNonuniformityNormalized: return this->GetRunLengthNonuniformityNormalized(); case LowGreyLevelRunEmphasis: return this->GetLowGreyLevelRunEmphasis(); case HighGreyLevelRunEmphasis: return this->GetHighGreyLevelRunEmphasis(); case ShortRunLowGreyLevelEmphasis: return this->GetShortRunLowGreyLevelEmphasis(); case ShortRunHighGreyLevelEmphasis: return this->GetShortRunHighGreyLevelEmphasis(); case LongRunLowGreyLevelEmphasis: return this->GetLongRunLowGreyLevelEmphasis(); case LongRunHighGreyLevelEmphasis: return this->GetLongRunHighGreyLevelEmphasis(); case RunPercentage: return this->GetRunPercentage(); case NumberOfRuns: return this->GetNumberOfRuns(); case GreyLevelVariance: return this->GetGreyLevelVariance(); case RunLengthVariance: return this->GetRunLengthVariance(); case RunEntropy: return this->GetRunEntropy(); default: return 0; } } template< typename THistogram> void EnhancedHistogramToRunLengthFeaturesFilter< THistogram>:: PrintSelf(std::ostream& os, Indent indent) const { Superclass::PrintSelf( os,indent ); } } // end of namespace Statistics } // end of namespace itk #endif diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp index fa8778bd9a..6b95762fda 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp @@ -1,604 +1,607 @@ #include // MITK #include #include #include // ITK #include #include #include #include // STL #include static void MatrixFeaturesTo(mitk::CoocurenceMatrixFeatures features, std::string prefix, mitk::GIFCooccurenceMatrix2::FeatureListType &featureList); static void CalculateMeanAndStdDevFeatures(std::vector featureList, mitk::CoocurenceMatrixFeatures &meanFeature, mitk::CoocurenceMatrixFeatures &stdFeature); static void NormalizeMatrixFeature(mitk::CoocurenceMatrixFeatures &features, std::size_t number); mitk::CoocurenceMatrixHolder::CoocurenceMatrixHolder(double min, double max, int number) : m_MinimumRange(min), m_MaximumRange(max), m_NumberOfBins(number) { m_Matrix.resize(number, number); m_Matrix.fill(0); m_Stepsize = (max - min) / (number); } int mitk::CoocurenceMatrixHolder::IntensityToIndex(double intensity) { return std::floor((intensity - m_MinimumRange) / m_Stepsize); } double mitk::CoocurenceMatrixHolder::IndexToMinIntensity(int index) { return m_MinimumRange + index * m_Stepsize; } double mitk::CoocurenceMatrixHolder::IndexToMeanIntensity(int index) { return m_MinimumRange + (index+0.5) * m_Stepsize; } double mitk::CoocurenceMatrixHolder::IndexToMaxIntensity(int index) { return m_MinimumRange + (index + 1) * m_Stepsize; } template void CalculateCoOcMatrix(itk::Image* itkImage, itk::Image* mask, itk::Offset offset, int range, mitk::CoocurenceMatrixHolder &holder) { typedef itk::Image ImageType; typedef itk::Image MaskImageType; typedef itk::ShapedNeighborhoodIterator ShapeIterType; typedef itk::ShapedNeighborhoodIterator ShapeMaskIterType; typedef itk::ImageRegionConstIterator ConstIterType; typedef itk::ImageRegionConstIterator ConstMaskIterType; itk::Size radius; radius.Fill(range+1); ShapeIterType imageOffsetIter(radius, itkImage, itkImage->GetLargestPossibleRegion()); ShapeMaskIterType maskOffsetIter(radius, mask, mask->GetLargestPossibleRegion()); imageOffsetIter.ActivateOffset(offset); maskOffsetIter.ActivateOffset(offset); ConstIterType imageIter(itkImage, itkImage->GetLargestPossibleRegion()); ConstMaskIterType maskIter(mask, mask->GetLargestPossibleRegion()); // iterator.GetIndex() + ci.GetNeighborhoodOffset() auto region = mask->GetLargestPossibleRegion(); while (!maskIter.IsAtEnd()) { auto ciMask = maskOffsetIter.Begin(); auto ciValue = imageOffsetIter.Begin(); if (maskIter.Value() > 0 && ciMask.Get() > 0 && imageIter.Get() == imageIter.Get() && ciValue.Get() == ciValue.Get() && region.IsInside(maskOffsetIter.GetIndex() + ciMask.GetNeighborhoodOffset())) { int i = holder.IntensityToIndex(imageIter.Get()); int j = holder.IntensityToIndex(ciValue.Get()); holder.m_Matrix(i, j) += 1; holder.m_Matrix(j, i) += 1; } ++imageOffsetIter; ++maskOffsetIter; ++imageIter; ++maskIter; } } void CalculateFeatures( mitk::CoocurenceMatrixHolder &holder, mitk::CoocurenceMatrixFeatures & results ) { auto pijMatrix = holder.m_Matrix; auto piMatrix = holder.m_Matrix; auto pjMatrix = holder.m_Matrix; double Ng = holder.m_NumberOfBins; int NgSize = holder.m_NumberOfBins; pijMatrix /= pijMatrix.sum(); piMatrix.rowwise().normalize(); pjMatrix.colwise().normalize(); for (int i = 0; i < holder.m_NumberOfBins; ++i) for (int j = 0; j < holder.m_NumberOfBins; ++j) { if (pijMatrix(i, j) != pijMatrix(i, j)) pijMatrix(i, j) = 0; if (piMatrix(i, j) != piMatrix(i, j)) piMatrix(i, j) = 0; if (pjMatrix(i, j) != pjMatrix(i, j)) pjMatrix(i, j) = 0; } Eigen::VectorXd piVector = pijMatrix.colwise().sum(); Eigen::VectorXd pjVector = pijMatrix.rowwise().sum(); double sigmai = 0;; for (int i = 0; i < holder.m_NumberOfBins; ++i) { - double iInt = holder.IndexToMeanIntensity(i); + double iInt = i + 1;// holder.IndexToMeanIntensity(i); results.RowAverage += iInt * piVector(i); if (piVector(i) > 0) { results.RowEntropy -= piVector(i) * std::log(piVector(i)) / std::log(2); } } for (int i = 0; i < holder.m_NumberOfBins; ++i) { - double iInt = holder.IndexToMeanIntensity(i); + double iInt = i + 1; // holder.IndexToMeanIntensity(i); results.RowVariance += (iInt - results.RowAverage)*(iInt - results.RowAverage) * piVector(i); } results.RowMaximum = piVector.maxCoeff(); sigmai = std::sqrt(results.RowVariance); Eigen::VectorXd pimj(NgSize); pimj.fill(0); Eigen::VectorXd pipj(2*NgSize); pipj.fill(0); results.JointMaximum += pijMatrix.maxCoeff(); for (int i = 0; i < holder.m_NumberOfBins; ++i) { for (int j = 0; j < holder.m_NumberOfBins; ++j) { - double iInt = holder.IndexToMeanIntensity(i); - double jInt = holder.IndexToMeanIntensity(j); + //double iInt = holder.IndexToMeanIntensity(i); + //double jInt = holder.IndexToMeanIntensity(j); + double iInt = i + 1;// holder.IndexToMeanIntensity(i); + double jInt = j + 1;// holder.IndexToMeanIntensity(j); double pij = pijMatrix(i, j); int deltaK = (i - j)>0?(i-j) : (j-i); pimj(deltaK) += pij; pipj(i + j) += pij; results.JointAverage += iInt * pij; if (pij > 0) { results.JointEntropy -= pij * std::log(pij) / std::log(2); results.FirstRowColumnEntropy -= pij * std::log(piVector(i)*pjVector(j)) / std::log(2); } if (piVector(i) > 0 && pjVector(j) > 0 ) { results.SecondRowColumnEntropy -= piVector(i)*pjVector(j) * std::log(piVector(i)*pjVector(j)) / std::log(2); } results.AngularSecondMoment += pij*pij; results.Contrast += (iInt - jInt)* (iInt - jInt) * pij; results.Dissimilarity += std::abs(iInt - jInt) * pij; results.InverseDifference += pij / (1 + (std::abs(iInt - jInt))); results.InverseDifferenceNormalised += pij / (1 + (std::abs(iInt - jInt) / Ng)); results.InverseDifferenceMoment += pij / (1 + (iInt - jInt)*(iInt - jInt)); results.InverseDifferenceMomentNormalised += pij / (1 + (iInt - jInt)*(iInt - jInt)/Ng/Ng); results.Autocorrelation += iInt*jInt * pij; double cluster = (iInt + jInt - 2 * results.RowAverage); results.ClusterTendency += cluster*cluster * pij; results.ClusterShade += cluster*cluster*cluster * pij; results.ClusterProminence += cluster*cluster*cluster*cluster * pij; if (iInt != jInt) { results.InverseVariance += pij / (iInt - jInt) / (iInt - jInt); } } } results.Correlation = 1 / sigmai / sigmai * (-results.RowAverage*results.RowAverage+ results.Autocorrelation); results.FirstMeasureOfInformationCorrelation = (results.JointEntropy - results.FirstRowColumnEntropy) / results.RowEntropy; if (results.JointEntropy < results.SecondRowColumnEntropy) { results.SecondMeasureOfInformationCorrelation = sqrt(1 - exp(-2 * (results.SecondRowColumnEntropy - results.JointEntropy))); } else { results.SecondMeasureOfInformationCorrelation = 0; } for (int i = 0; i < holder.m_NumberOfBins; ++i) { for (int j = 0; j < holder.m_NumberOfBins; ++j) { - double iInt = holder.IndexToMeanIntensity(i); + //double iInt = holder.IndexToMeanIntensity(i); //double jInt = holder.IndexToMeanIntensity(j); + double iInt = i + 1; double pij = pijMatrix(i, j); results.JointVariance += (iInt - results.JointAverage)* (iInt - results.JointAverage)*pij; } } for (int k = 0; k < NgSize; ++k) { results.DifferenceAverage += k* pimj(k); if (pimj(k) > 0) { results.DifferenceEntropy -= pimj(k) * log(pimj(k)) / std::log(2); } } for (int k = 0; k < NgSize; ++k) { results.DifferenceVariance += (results.DifferenceAverage-k)* (results.DifferenceAverage-k)*pimj(k); } for (int k = 0; k <2* NgSize ; ++k) { results.SumAverage += (2+k)* pipj(k); if (pipj(k) > 0) { results.SumEntropy -= pipj(k) * log(pipj(k)) / std::log(2); } } for (int k = 0; k < 2*NgSize; ++k) { results.SumVariance += (2+k - results.SumAverage)* (2+k - results.SumAverage)*pipj(k); } //MITK_INFO << std::endl << holder.m_Matrix; //MITK_INFO << std::endl << pijMatrix; //MITK_INFO << std::endl << piMatrix; //MITK_INFO << std::endl << pjMatrix; //for (int i = 0; i < holder.m_NumberOfBins; ++i) //{ // MITK_INFO << "Bin " << i << " Min: " << holder.IndexToMinIntensity(i) << " Max: " << holder.IndexToMaxIntensity(i); //} //MITK_INFO << pimj; //MITK_INFO << pipj; } template void CalculateCoocurenceFeatures(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFCooccurenceMatrix2::FeatureListType & featureList, mitk::GIFCooccurenceMatrix2::GIFCooccurenceMatrix2Configuration config) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::MinimumMaximumImageCalculator MinMaxComputerType; typedef itk::Neighborhood NeighborhoodType; typedef itk::Offset OffsetType; /////////////////////////////////////////////////////////////////////////////////////////////// typename MinMaxComputerType::Pointer minMaxComputer = MinMaxComputerType::New(); minMaxComputer->SetImage(itkImage); minMaxComputer->Compute(); double rangeMin = -0.5 + minMaxComputer->GetMinimum(); double rangeMax = 0.5 + minMaxComputer->GetMaximum(); if (config.UseMinimumIntensity) rangeMin = config.MinimumIntensity; if (config.UseMaximumIntensity) rangeMax = config.MaximumIntensity; // Define Range int numberOfBins = config.Bins; typename MaskType::Pointer maskImage = MaskType::New(); mitk::CastToItkImage(mask, maskImage); //Find possible directions std::vector < itk::Offset > offsetVector; NeighborhoodType hood; hood.SetRadius(1); unsigned int centerIndex = hood.GetCenterNeighborhoodIndex(); OffsetType offset; for (unsigned int d = 0; d < centerIndex; d++) { offset = hood.GetOffset(d); bool useOffset = true; for (unsigned int i = 0; i < VImageDimension; ++i) { offset[i] *= config.range; if (config.direction == i + 2 && offset[i] != 0) { useOffset = false; } } if (useOffset) { offsetVector.push_back(offset); } } if (config.direction == 1) { offsetVector.clear(); offset[0] = 0; offset[1] = 0; offset[2] = 1; } std::vector resultVector; mitk::CoocurenceMatrixHolder holderOverall(rangeMin, rangeMax, numberOfBins); mitk::CoocurenceMatrixFeatures overallFeature; for (std::size_t i = 0; i < offsetVector.size(); ++i) { if (config.direction > 1) { if (offsetVector[i][config.direction - 2] != 0) { continue; } } offset = offsetVector[i]; mitk::CoocurenceMatrixHolder holder(rangeMin, rangeMax, numberOfBins); mitk::CoocurenceMatrixFeatures coocResults; CalculateCoOcMatrix(itkImage, maskImage, offset, config.range, holder); holderOverall.m_Matrix += holder.m_Matrix; CalculateFeatures(holder, coocResults); resultVector.push_back(coocResults); } CalculateFeatures(holderOverall, overallFeature); //NormalizeMatrixFeature(overallFeature, offsetVector.size()); mitk::CoocurenceMatrixFeatures featureMean; mitk::CoocurenceMatrixFeatures featureStd; CalculateMeanAndStdDevFeatures(resultVector, featureMean, featureStd); std::ostringstream ss; ss << config.range; std::string strRange = ss.str(); MatrixFeaturesTo(overallFeature, "co-occ. 2 (" + strRange + ") overall", featureList); MatrixFeaturesTo(featureMean, "co-occ. 2 (" + strRange + ") mean", featureList); MatrixFeaturesTo(featureStd, "co-occ. 2 (" + strRange + ") std.dev.", featureList); } static void MatrixFeaturesTo(mitk::CoocurenceMatrixFeatures features, std::string prefix, mitk::GIFCooccurenceMatrix2::FeatureListType &featureList) { featureList.push_back(std::make_pair(prefix + " Joint Maximum", features.JointMaximum)); featureList.push_back(std::make_pair(prefix + " Joint Average", features.JointAverage)); featureList.push_back(std::make_pair(prefix + " Joint Variance", features.JointVariance)); featureList.push_back(std::make_pair(prefix + " Joint Entropy", features.JointEntropy)); featureList.push_back(std::make_pair(prefix + " Row Maximum", features.RowMaximum)); featureList.push_back(std::make_pair(prefix + " Row Average", features.RowAverage)); featureList.push_back(std::make_pair(prefix + " Row Variance", features.RowVariance)); featureList.push_back(std::make_pair(prefix + " Row Entropy", features.RowEntropy)); featureList.push_back(std::make_pair(prefix + " First Row-Column Entropy", features.FirstRowColumnEntropy)); featureList.push_back(std::make_pair(prefix + " Second Row-Column Entropy", features.SecondRowColumnEntropy)); featureList.push_back(std::make_pair(prefix + " Difference Average", features.DifferenceAverage)); featureList.push_back(std::make_pair(prefix + " Difference Variance", features.DifferenceVariance)); featureList.push_back(std::make_pair(prefix + " Difference Entropy", features.DifferenceEntropy)); featureList.push_back(std::make_pair(prefix + " Sum Average", features.SumAverage)); featureList.push_back(std::make_pair(prefix + " Sum Variance", features.SumVariance)); featureList.push_back(std::make_pair(prefix + " Sum Entropy", features.SumEntropy)); featureList.push_back(std::make_pair(prefix + " Angular Second Moment", features.AngularSecondMoment)); featureList.push_back(std::make_pair(prefix + " Contrast", features.Contrast)); featureList.push_back(std::make_pair(prefix + " Dissimilarity", features.Dissimilarity)); featureList.push_back(std::make_pair(prefix + " Inverse Difference", features.InverseDifference)); featureList.push_back(std::make_pair(prefix + " Inverse Difference Normalized", features.InverseDifferenceNormalised)); featureList.push_back(std::make_pair(prefix + " Inverse Difference Moment", features.InverseDifferenceMoment)); featureList.push_back(std::make_pair(prefix + " Inverse Difference Moment Normalized", features.InverseDifferenceMomentNormalised)); featureList.push_back(std::make_pair(prefix + " Inverse Variance", features.InverseVariance)); featureList.push_back(std::make_pair(prefix + " Correlation", features.Correlation)); featureList.push_back(std::make_pair(prefix + " Autocorrleation", features.Autocorrelation)); featureList.push_back(std::make_pair(prefix + " Cluster Tendency", features.ClusterTendency)); featureList.push_back(std::make_pair(prefix + " Cluster Shade", features.ClusterShade)); featureList.push_back(std::make_pair(prefix + " Cluster Prominence", features.ClusterProminence)); featureList.push_back(std::make_pair(prefix + " First Measure of Information Correlation", features.FirstMeasureOfInformationCorrelation)); featureList.push_back(std::make_pair(prefix + " Second Measure of Information Correlation", features.SecondMeasureOfInformationCorrelation)); } static void CalculateMeanAndStdDevFeatures(std::vector featureList, mitk::CoocurenceMatrixFeatures &meanFeature, mitk::CoocurenceMatrixFeatures &stdFeature) { #define ADDFEATURE(a) meanFeature.a += featureList[i].a;stdFeature.a += featureList[i].a*featureList[i].a #define CALCVARIANCE(a) stdFeature.a =sqrt(stdFeature.a - meanFeature.a*meanFeature.a) for (std::size_t i = 0; i < featureList.size(); ++i) { ADDFEATURE(JointMaximum); ADDFEATURE(JointAverage); ADDFEATURE(JointVariance); ADDFEATURE(JointEntropy); ADDFEATURE(RowMaximum); ADDFEATURE(RowAverage); ADDFEATURE(RowVariance); ADDFEATURE(RowEntropy); ADDFEATURE(FirstRowColumnEntropy); ADDFEATURE(SecondRowColumnEntropy); ADDFEATURE(DifferenceAverage); ADDFEATURE(DifferenceVariance); ADDFEATURE(DifferenceEntropy); ADDFEATURE(SumAverage); ADDFEATURE(SumVariance); ADDFEATURE(SumEntropy); ADDFEATURE(AngularSecondMoment); ADDFEATURE(Contrast); ADDFEATURE(Dissimilarity); ADDFEATURE(InverseDifference); ADDFEATURE(InverseDifferenceNormalised); ADDFEATURE(InverseDifferenceMoment); ADDFEATURE(InverseDifferenceMomentNormalised); ADDFEATURE(InverseVariance); ADDFEATURE(Correlation); ADDFEATURE(Autocorrelation); ADDFEATURE(ClusterShade); ADDFEATURE(ClusterTendency); ADDFEATURE(ClusterProminence); ADDFEATURE(FirstMeasureOfInformationCorrelation); ADDFEATURE(SecondMeasureOfInformationCorrelation); } NormalizeMatrixFeature(meanFeature, featureList.size()); NormalizeMatrixFeature(stdFeature, featureList.size()); CALCVARIANCE(JointMaximum); CALCVARIANCE(JointAverage); CALCVARIANCE(JointVariance); CALCVARIANCE(JointEntropy); CALCVARIANCE(RowMaximum); CALCVARIANCE(RowAverage); CALCVARIANCE(RowVariance); CALCVARIANCE(RowEntropy); CALCVARIANCE(FirstRowColumnEntropy); CALCVARIANCE(SecondRowColumnEntropy); CALCVARIANCE(DifferenceAverage); CALCVARIANCE(DifferenceVariance); CALCVARIANCE(DifferenceEntropy); CALCVARIANCE(SumAverage); CALCVARIANCE(SumVariance); CALCVARIANCE(SumEntropy); CALCVARIANCE(AngularSecondMoment); CALCVARIANCE(Contrast); CALCVARIANCE(Dissimilarity); CALCVARIANCE(InverseDifference); CALCVARIANCE(InverseDifferenceNormalised); CALCVARIANCE(InverseDifferenceMoment); CALCVARIANCE(InverseDifferenceMomentNormalised); CALCVARIANCE(InverseVariance); CALCVARIANCE(Correlation); CALCVARIANCE(Autocorrelation); CALCVARIANCE(ClusterShade); CALCVARIANCE(ClusterTendency); CALCVARIANCE(ClusterProminence); CALCVARIANCE(FirstMeasureOfInformationCorrelation); CALCVARIANCE(SecondMeasureOfInformationCorrelation); #undef ADDFEATURE #undef CALCVARIANCE } static void NormalizeMatrixFeature(mitk::CoocurenceMatrixFeatures &features, std::size_t number) { features.JointMaximum = features.JointMaximum / number; features.JointAverage = features.JointAverage / number; features.JointVariance = features.JointVariance / number; features.JointEntropy = features.JointEntropy / number; features.RowMaximum = features.RowMaximum / number; features.RowAverage = features.RowAverage / number; features.RowVariance = features.RowVariance / number; features.RowEntropy = features.RowEntropy / number; features.FirstRowColumnEntropy = features.FirstRowColumnEntropy / number; features.SecondRowColumnEntropy = features.SecondRowColumnEntropy / number; features.DifferenceAverage = features.DifferenceAverage / number; features.DifferenceVariance = features.DifferenceVariance / number; features.DifferenceEntropy = features.DifferenceEntropy / number; features.SumAverage = features.SumAverage / number; features.SumVariance = features.SumVariance / number; features.SumEntropy = features.SumEntropy / number; features.AngularSecondMoment = features.AngularSecondMoment / number; features.Contrast = features.Contrast / number; features.Dissimilarity = features.Dissimilarity / number; features.InverseDifference = features.InverseDifference / number; features.InverseDifferenceNormalised = features.InverseDifferenceNormalised / number; features.InverseDifferenceMoment = features.InverseDifferenceMoment / number; features.InverseDifferenceMomentNormalised = features.InverseDifferenceMomentNormalised / number; features.InverseVariance = features.InverseVariance / number; features.Correlation = features.Correlation / number; features.Autocorrelation = features.Autocorrelation / number; features.ClusterShade = features.ClusterShade / number; features.ClusterTendency = features.ClusterTendency / number; features.ClusterProminence = features.ClusterProminence / number; features.FirstMeasureOfInformationCorrelation = features.FirstMeasureOfInformationCorrelation / number; features.SecondMeasureOfInformationCorrelation = features.SecondMeasureOfInformationCorrelation / number; } mitk::GIFCooccurenceMatrix2::GIFCooccurenceMatrix2(): m_Range(1.0), m_Bins(128) { SetShortName("cooc2"); SetLongName("cooccurence2"); } mitk::GIFCooccurenceMatrix2::FeatureListType mitk::GIFCooccurenceMatrix2::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; GIFCooccurenceMatrix2Configuration config; config.direction = GetDirection(); config.range = m_Range; config.MinimumIntensity = GetMinimumIntensity(); config.MaximumIntensity = GetMaximumIntensity(); config.UseMinimumIntensity = GetUseMinimumIntensity(); config.UseMaximumIntensity = GetUseMaximumIntensity(); config.Bins = GetBins(); AccessByItk_3(image, CalculateCoocurenceFeatures, mask, featureList,config); return featureList; } mitk::GIFCooccurenceMatrix2::FeatureNameListType mitk::GIFCooccurenceMatrix2::GetFeatureNames() { FeatureNameListType featureList; return featureList; } void mitk::GIFCooccurenceMatrix2::AddArguments(mitkCommandLineParser &parser) { std::string name = GetOptionPrefix(); parser.addArgument(GetLongName(), name, mitkCommandLineParser::String, "Use Co-occurence matrix", "calculates Co-occurence based features (new implementation)", us::Any()); parser.addArgument(name+"::range", name+"::range", mitkCommandLineParser::String, "Cooc 2 Range", "Define the range that is used (Semicolon-separated)", us::Any()); parser.addArgument(name + "::bins", name + "::bins", mitkCommandLineParser::String, "Cooc 2 Number of Bins", "Define the number of bins that is used ", us::Any()); } void mitk::GIFCooccurenceMatrix2::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &, const Image::Pointer &maskNoNAN, FeatureListType &featureList) { auto parsedArgs = GetParameter(); std::string name = GetOptionPrefix(); if (parsedArgs.count(GetLongName())) { std::vector ranges; if (parsedArgs.count(name + "::range")) { ranges = SplitDouble(parsedArgs[name + "::range"].ToString(), ';'); } else { ranges.push_back(1); } if (parsedArgs.count(name + "::bins")) { auto bins = SplitDouble(parsedArgs[name + "::bins"].ToString(), ';')[0]; this->SetBins(bins); } for (std::size_t i = 0; i < ranges.size(); ++i) { MITK_INFO << "Start calculating coocurence with range " << ranges[i] << "...."; this->SetRange(ranges[i]); auto localResults = this->CalculateFeatures(feature, maskNoNAN); featureList.insert(featureList.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating coocurence with range " << ranges[i] << "...."; } } } diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderStatistics.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderStatistics.cpp index 114656a8fa..119fefd241 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderStatistics.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderStatistics.cpp @@ -1,301 +1,314 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include // MITK #include #include #include // ITK #include #include // STL #include template void CalculateFirstOrderStatistics(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFFirstOrderStatistics::FeatureListType & featureList, mitk::GIFFirstOrderStatistics::ParameterStruct params) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::LabelStatisticsImageFilter FilterType; typedef typename FilterType::HistogramType HistogramType; typedef typename HistogramType::IndexType HIndexType; typedef itk::MinimumMaximumImageCalculator MinMaxComputerType; typename MaskType::Pointer maskImage = MaskType::New(); mitk::CastToItkImage(mask, maskImage); typename MinMaxComputerType::Pointer minMaxComputer = MinMaxComputerType::New(); minMaxComputer->SetImage(itkImage); minMaxComputer->Compute(); double imageRange = minMaxComputer->GetMaximum() - minMaxComputer->GetMinimum(); typename FilterType::Pointer labelStatisticsImageFilter = FilterType::New(); labelStatisticsImageFilter->SetInput( itkImage ); labelStatisticsImageFilter->SetLabelInput(maskImage); labelStatisticsImageFilter->SetUseHistograms(true); if (params.m_UseCtRange) { labelStatisticsImageFilter->SetHistogramParameters(1024.5+3096.5, -1024.5,3096.5); } else { double min = minMaxComputer->GetMinimum() -0.5; double max = minMaxComputer->GetMaximum() +0.5; if (params.UseMinimumIntensity) min = params.MinimumIntensity; if (params.UseMaximumIntensity) max = params.MaximumIntensity; labelStatisticsImageFilter->SetHistogramParameters(params.Bins, min,max); } labelStatisticsImageFilter->Update(); // --------------- Range -------------------- double range = labelStatisticsImageFilter->GetMaximum(1) - labelStatisticsImageFilter->GetMinimum(1); // --------------- Uniformity, Entropy -------------------- double count = labelStatisticsImageFilter->GetCount(1); //double std_dev = labelStatisticsImageFilter->GetSigma(1); - double uncorrected_std_dev = std::sqrt((count - 1) / count * labelStatisticsImageFilter->GetVariance(1)); double mean = labelStatisticsImageFilter->GetMean(1); auto histogram = labelStatisticsImageFilter->GetHistogram(1); bool histogramIsCalculated = histogram; HIndexType index; index.SetSize(1); double uniformity = 0; double entropy = 0; double squared_sum = 0; double kurtosis = 0; double mean_absolut_deviation = 0; double skewness = 0; double sum_prob = 0; double binWidth = 0; double p05th = 0, p10th = 0, p15th = 0, p20th = 0, p25th = 0, p30th = 0, p35th = 0, p40th = 0, p45th = 0, p50th = 0; double p55th = 0, p60th = 0, p65th = 0, p70th = 0, p75th = 0, p80th = 0, p85th = 0, p90th = 0, p95th = 0; double voxelValue = 0; if (histogramIsCalculated) { binWidth = histogram->GetBinMax(0, 0) - histogram->GetBinMin(0, 0); p05th = histogram->Quantile(0, 0.05); p10th = histogram->Quantile(0, 0.10); p15th = histogram->Quantile(0, 0.15); p20th = histogram->Quantile(0, 0.20); p25th = histogram->Quantile(0, 0.25); p30th = histogram->Quantile(0, 0.30); p35th = histogram->Quantile(0, 0.35); p40th = histogram->Quantile(0, 0.40); p45th = histogram->Quantile(0, 0.45); p50th = histogram->Quantile(0, 0.50); p55th = histogram->Quantile(0, 0.55); p60th = histogram->Quantile(0, 0.60); p65th = histogram->Quantile(0, 0.65); p70th = histogram->Quantile(0, 0.70); p75th = histogram->Quantile(0, 0.75); p80th = histogram->Quantile(0, 0.80); p85th = histogram->Quantile(0, 0.85); p90th = histogram->Quantile(0, 0.90); p95th = histogram->Quantile(0, 0.95); } double Log2=log(2); + double mode_bin; + double mode_value = 0; + double variance = 0; if (histogramIsCalculated) { for (int i = 0; i < (int)(histogram->GetSize(0)); ++i) { index[0] = i; double prob = histogram->GetFrequency(index); - if (prob < 0.1) + + if (prob < 0.00000001) continue; - if (histogramIsCalculated) + + voxelValue = histogram->GetBinMin(0, i) + binWidth * 0.5; + + if (prob > mode_value) { - voxelValue = histogram->GetBinMin(0, i) + binWidth * 0.5; + mode_value = prob; + mode_bin = voxelValue; } + sum_prob += prob; squared_sum += prob * voxelValue*voxelValue; prob /= count; mean_absolut_deviation += prob* std::abs(voxelValue - mean); + variance += prob * (voxelValue - mean) * (voxelValue - mean); kurtosis += prob* (voxelValue - mean) * (voxelValue - mean) * (voxelValue - mean) * (voxelValue - mean); skewness += prob* (voxelValue - mean) * (voxelValue - mean) * (voxelValue - mean); uniformity += prob*prob; if (prob > 0) { entropy += prob * std::log(prob) / Log2; } } } entropy = -entropy; + double uncorrected_std_dev = std::sqrt(variance); double rms = std::sqrt(squared_sum / count); - kurtosis = kurtosis / (uncorrected_std_dev*uncorrected_std_dev * uncorrected_std_dev*uncorrected_std_dev); - skewness = skewness / (uncorrected_std_dev*uncorrected_std_dev * uncorrected_std_dev); + kurtosis = kurtosis / (variance * variance); + skewness = skewness / (variance * uncorrected_std_dev); //mean_absolut_deviation = mean_absolut_deviation; double coveredGrayValueRange = range / imageRange; featureList.push_back(std::make_pair("FirstOrder Range",range)); featureList.push_back(std::make_pair("FirstOrder Uniformity",uniformity)); featureList.push_back(std::make_pair("FirstOrder Entropy",entropy)); featureList.push_back(std::make_pair("FirstOrder Energy",squared_sum)); featureList.push_back(std::make_pair("FirstOrder RMS",rms)); featureList.push_back(std::make_pair("FirstOrder Kurtosis", kurtosis)); featureList.push_back(std::make_pair("FirstOrder Excess Kurtosis", kurtosis-3)); featureList.push_back(std::make_pair("FirstOrder Skewness",skewness)); featureList.push_back(std::make_pair("FirstOrder Mean absolute deviation",mean_absolut_deviation)); featureList.push_back(std::make_pair("FirstOrder Covered Image Intensity Range",coveredGrayValueRange)); featureList.push_back(std::make_pair("FirstOrder Minimum",labelStatisticsImageFilter->GetMinimum(1))); featureList.push_back(std::make_pair("FirstOrder Maximum",labelStatisticsImageFilter->GetMaximum(1))); featureList.push_back(std::make_pair("FirstOrder Mean",labelStatisticsImageFilter->GetMean(1))); - featureList.push_back(std::make_pair("FirstOrder Variance",labelStatisticsImageFilter->GetVariance(1))); + featureList.push_back(std::make_pair("FirstOrder Unbiased Variance", labelStatisticsImageFilter->GetVariance(1))); //Siehe Definition von Unbiased Variance estimation. (Wird nicht durch n sondern durch n-1 normalisiert) + featureList.push_back(std::make_pair("FirstOrder Biased Variance", variance)); featureList.push_back(std::make_pair("FirstOrder Sum",labelStatisticsImageFilter->GetSum(1))); - featureList.push_back(std::make_pair("FirstOrder Median",labelStatisticsImageFilter->GetMedian(1))); + featureList.push_back(std::make_pair("FirstOrder Median", labelStatisticsImageFilter->GetMedian(1))); + featureList.push_back(std::make_pair("FirstOrder Mode", mode_bin)); + featureList.push_back(std::make_pair("FirstOrder Mode Probability", mode_value)); featureList.push_back(std::make_pair("FirstOrder Standard deviation",labelStatisticsImageFilter->GetSigma(1))); featureList.push_back(std::make_pair("FirstOrder No. of Voxel",labelStatisticsImageFilter->GetCount(1))); featureList.push_back(std::make_pair("FirstOrder 05th Percentile", p05th)); featureList.push_back(std::make_pair("FirstOrder 10th Percentile", p10th)); featureList.push_back(std::make_pair("FirstOrder 15th Percentile", p15th)); featureList.push_back(std::make_pair("FirstOrder 20th Percentile", p20th)); featureList.push_back(std::make_pair("FirstOrder 25th Percentile", p25th)); featureList.push_back(std::make_pair("FirstOrder 30th Percentile", p30th)); featureList.push_back(std::make_pair("FirstOrder 35th Percentile", p35th)); featureList.push_back(std::make_pair("FirstOrder 40th Percentile", p40th)); featureList.push_back(std::make_pair("FirstOrder 45th Percentile", p45th)); featureList.push_back(std::make_pair("FirstOrder 50th Percentile", p50th)); featureList.push_back(std::make_pair("FirstOrder 55th Percentile", p55th)); featureList.push_back(std::make_pair("FirstOrder 60th Percentile", p60th)); featureList.push_back(std::make_pair("FirstOrder 65th Percentile", p65th)); featureList.push_back(std::make_pair("FirstOrder 70th Percentile", p70th)); featureList.push_back(std::make_pair("FirstOrder 75th Percentile", p75th)); featureList.push_back(std::make_pair("FirstOrder 80th Percentile", p80th)); featureList.push_back(std::make_pair("FirstOrder 85th Percentile", p85th)); featureList.push_back(std::make_pair("FirstOrder 90th Percentile", p90th)); featureList.push_back(std::make_pair("FirstOrder 95th Percentile", p95th)); featureList.push_back(std::make_pair("FirstOrder Interquartile Range",(p75th - p25th))); //Calculate the robust mean absolute deviation //First, set all frequencies to 0 that are <10th or >90th percentile double meanRobust = 0.0; double robustMeanAbsoluteDeviation = 0.0; if (histogramIsCalculated) { for (int i = 0; i < (int)(histogram->GetSize(0)); ++i) { index[0] = i; if (histogram->GetBinMax(0, i) < p10th) { histogram->SetFrequencyOfIndex(index, 0); } else if (histogram->GetBinMin(0, i) > p90th) { histogram->SetFrequencyOfIndex(index, 0); } } //Calculate the mean for (int i = 0; i < (int)(histogram->GetSize(0)); ++i) { index[0] = i; meanRobust += histogram->GetFrequency(index) * 0.5 * (histogram->GetBinMin(0, i) + histogram->GetBinMax(0, i)); } meanRobust = meanRobust / histogram->GetTotalFrequency(); for (int i = 0; i < (int)(histogram->GetSize(0)); ++i) { index[0] = i; robustMeanAbsoluteDeviation += std::abs(histogram->GetFrequency(index) * ((0.5 * (histogram->GetBinMin(0, i) + histogram->GetBinMax(0, i))) - meanRobust )); } robustMeanAbsoluteDeviation = robustMeanAbsoluteDeviation / histogram->GetTotalFrequency(); } featureList.push_back(std::make_pair("FirstOrder Robust Mean", meanRobust)); featureList.push_back(std::make_pair("FirstOrder Robust Mean Absolute Deviation", robustMeanAbsoluteDeviation)); } mitk::GIFFirstOrderStatistics::GIFFirstOrderStatistics() : m_HistogramSize(256), m_UseCtRange(false) { SetShortName("fo"); SetLongName("first-order"); } mitk::GIFFirstOrderStatistics::FeatureListType mitk::GIFFirstOrderStatistics::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; ParameterStruct params; params.m_HistogramSize = this->m_HistogramSize; params.m_UseCtRange = this->m_UseCtRange; params.MinimumIntensity = GetMinimumIntensity(); params.MaximumIntensity = GetMaximumIntensity(); params.UseMinimumIntensity = GetUseMinimumIntensity(); params.UseMaximumIntensity = GetUseMaximumIntensity(); params.Bins = GetBins(); AccessByItk_3(image, CalculateFirstOrderStatistics, mask, featureList, params); return featureList; } mitk::GIFFirstOrderStatistics::FeatureNameListType mitk::GIFFirstOrderStatistics::GetFeatureNames() { FeatureNameListType featureList; featureList.push_back("FirstOrder Minimum"); featureList.push_back("FirstOrder Maximum"); featureList.push_back("FirstOrder Mean"); featureList.push_back("FirstOrder Variance"); featureList.push_back("FirstOrder Sum"); featureList.push_back("FirstOrder Median"); featureList.push_back("FirstOrder Standard deviation"); featureList.push_back("FirstOrder No. of Voxel"); return featureList; } void mitk::GIFFirstOrderStatistics::AddArguments(mitkCommandLineParser &parser) { std::string name = GetOptionPrefix(); parser.addArgument(GetLongName(), name, mitkCommandLineParser::String, "Use Volume-Statistic", "calculates volume based features", us::Any()); } void mitk::GIFFirstOrderStatistics::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &, const Image::Pointer &maskNoNAN, FeatureListType &featureList) { auto parsedArgs = GetParameter(); if (parsedArgs.count(GetLongName())) { MITK_INFO << "Start calculating first order features ...."; auto localResults = this->CalculateFeatures(feature, maskNoNAN); featureList.insert(featureList.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating first order features...."; } } diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelRunLength.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelRunLength.cpp index 219a91f8bf..eb446f7f21 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelRunLength.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelRunLength.cpp @@ -1,341 +1,340 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include // MITK #include #include #include // ITK #include #include // STL #include template void CalculateGrayLevelRunLengthFeatures(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFGrayLevelRunLength::FeatureListType & featureList, mitk::GIFGrayLevelRunLength::ParameterStruct params) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::Statistics::EnhancedScalarImageToRunLengthFeaturesFilter FilterType; typedef itk::MinimumMaximumImageCalculator MinMaxComputerType; typedef typename FilterType::RunLengthFeaturesFilterType TextureFilterType; typename MaskType::Pointer maskImage = MaskType::New(); mitk::CastToItkImage(mask, maskImage); typename FilterType::Pointer filter = FilterType::New(); typename FilterType::Pointer filter2 = FilterType::New(); typename FilterType::OffsetVector::Pointer newOffset = FilterType::OffsetVector::New(); auto oldOffsets = filter->GetOffsets(); auto oldOffsetsIterator = oldOffsets->Begin(); while (oldOffsetsIterator != oldOffsets->End()) { bool continueOuterLoop = false; typename FilterType::OffsetType offset = oldOffsetsIterator->Value(); for (unsigned int i = 0; i < VImageDimension; ++i) { if (params.m_Direction == i + 2 && offset[i] != 0) { continueOuterLoop = true; } } if (params.m_Direction == 1) { offset[0] = 0; offset[1] = 0; offset[2] = 1; newOffset->push_back(offset); break; } oldOffsetsIterator++; if (continueOuterLoop) continue; newOffset->push_back(offset); } filter->SetOffsets(newOffset); filter2->SetOffsets(newOffset); // All features are required typename FilterType::FeatureNameVectorPointer requestedFeatures = FilterType::FeatureNameVector::New(); requestedFeatures->push_back(TextureFilterType::ShortRunEmphasis); requestedFeatures->push_back(TextureFilterType::LongRunEmphasis); requestedFeatures->push_back(TextureFilterType::GreyLevelNonuniformity); requestedFeatures->push_back(TextureFilterType::GreyLevelNonuniformityNormalized); requestedFeatures->push_back(TextureFilterType::RunLengthNonuniformity); requestedFeatures->push_back(TextureFilterType::RunLengthNonuniformityNormalized); requestedFeatures->push_back(TextureFilterType::LowGreyLevelRunEmphasis); requestedFeatures->push_back(TextureFilterType::HighGreyLevelRunEmphasis); requestedFeatures->push_back(TextureFilterType::ShortRunLowGreyLevelEmphasis); requestedFeatures->push_back(TextureFilterType::ShortRunHighGreyLevelEmphasis); requestedFeatures->push_back(TextureFilterType::LongRunLowGreyLevelEmphasis); requestedFeatures->push_back(TextureFilterType::LongRunHighGreyLevelEmphasis); requestedFeatures->push_back(TextureFilterType::RunPercentage); requestedFeatures->push_back(TextureFilterType::NumberOfRuns); requestedFeatures->push_back(TextureFilterType::GreyLevelVariance); requestedFeatures->push_back(TextureFilterType::RunLengthVariance); requestedFeatures->push_back(TextureFilterType::RunEntropy); typename MinMaxComputerType::Pointer minMaxComputer = MinMaxComputerType::New(); minMaxComputer->SetImage(itkImage); minMaxComputer->Compute(); filter->SetInput(itkImage); filter->SetMaskImage(maskImage); filter->SetRequestedFeatures(requestedFeatures); filter2->SetInput(itkImage); filter2->SetMaskImage(maskImage); filter2->SetRequestedFeatures(requestedFeatures); int rangeOfPixels = params.Bins; if (rangeOfPixels < 2) rangeOfPixels = 256; if (params.m_UseCtRange) { filter->SetPixelValueMinMax((TPixel)(-1024.5),(TPixel)(3096.5)); filter->SetNumberOfBinsPerAxis(3096.5 + 1024.5); filter2->SetPixelValueMinMax((TPixel)(-1024.5), (TPixel)(3096.5)); filter2->SetNumberOfBinsPerAxis(3096.5 + 1024.5); } else { double minRange = minMaxComputer->GetMinimum() - 0.5; double maxRange = minMaxComputer->GetMaximum() + 0.5; + if (params.UseMinimumIntensity) minRange = params.MinimumIntensity; if (params.UseMaximumIntensity) maxRange = params.MaximumIntensity; + MITK_INFO << minRange << " " << maxRange; + filter->SetPixelValueMinMax(minRange, maxRange); filter->SetNumberOfBinsPerAxis(rangeOfPixels); filter2->SetPixelValueMinMax(minRange, maxRange); filter2->SetNumberOfBinsPerAxis(rangeOfPixels); } filter->SetDistanceValueMinMax(0, rangeOfPixels); filter2->SetDistanceValueMinMax(0, rangeOfPixels); filter2->CombinedFeatureCalculationOn(); filter->Update(); filter2->Update(); auto featureMeans = filter->GetFeatureMeans (); auto featureStd = filter->GetFeatureStandardDeviations(); auto featureCombined = filter2->GetFeatureMeans(); std::ostringstream ss; ss << rangeOfPixels; std::string strRange = ss.str(); for (std::size_t i = 0; i < featureMeans->size(); ++i) { switch (i) { case TextureFilterType::ShortRunEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") ShortRunEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") ShortRunEmphasis Std.", featureStd->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") ShortRunEmphasis Comb.", featureCombined->ElementAt(i))); break; case TextureFilterType::LongRunEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") LongRunEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") LongRunEmphasis Std.", featureStd->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") LongRunEmphasis Comb.", featureCombined->ElementAt(i))); break; case TextureFilterType::GreyLevelNonuniformity : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") GreyLevelNonuniformity Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") GreyLevelNonuniformity Std.", featureStd->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") GreyLevelNonuniformity Comb.", featureCombined->ElementAt(i))); break; case TextureFilterType::GreyLevelNonuniformityNormalized : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") GreyLevelNonuniformityNormalized Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") GreyLevelNonuniformityNormalized Std.", featureStd->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") GreyLevelNonuniformityNormalized Comb.", featureCombined->ElementAt(i))); break; case TextureFilterType::RunLengthNonuniformity : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunLengthNonuniformity Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") RunLengthNonuniformity Std.", featureStd->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") RunLengthNonuniformity Comb.", featureCombined->ElementAt(i))); break; case TextureFilterType::RunLengthNonuniformityNormalized : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunLengthNonuniformityNormalized Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") RunLengthNonuniformityNormalized Std.", featureStd->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") RunLengthNonuniformityNormalized Comb.", featureCombined->ElementAt(i))); break; case TextureFilterType::LowGreyLevelRunEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") LowGreyLevelRunEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") LowGreyLevelRunEmphasis Std.", featureStd->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") LowGreyLevelRunEmphasis Comb.", featureCombined->ElementAt(i))); break; case TextureFilterType::HighGreyLevelRunEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") HighGreyLevelRunEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") HighGreyLevelRunEmphasis Std.", featureStd->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") HighGreyLevelRunEmphasis Comb.", featureCombined->ElementAt(i))); break; case TextureFilterType::ShortRunLowGreyLevelEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") ShortRunLowGreyLevelEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") ShortRunLowGreyLevelEmphasis Std.", featureStd->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") ShortRunLowGreyLevelEmphasis Comb.", featureCombined->ElementAt(i))); break; case TextureFilterType::ShortRunHighGreyLevelEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") ShortRunHighGreyLevelEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") ShortRunHighGreyLevelEmphasis Std.", featureStd->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") ShortRunHighGreyLevelEmphasis Comb.", featureCombined->ElementAt(i))); break; case TextureFilterType::LongRunLowGreyLevelEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") LongRunLowGreyLevelEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") LongRunLowGreyLevelEmphasis Std.", featureStd->ElementAt(i))); - featureList.push_back(std::make_pair("RunLength. (" + strRange + ") LongRunLowGreyLevelEmphasis Comb.", featureStd->ElementAt(i))); + featureList.push_back(std::make_pair("RunLength. (" + strRange + ") LongRunLowGreyLevelEmphasis Comb.", featureCombined->ElementAt(i))); break; case TextureFilterType::LongRunHighGreyLevelEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") LongRunHighGreyLevelEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") LongRunHighGreyLevelEmphasis Std.", featureStd->ElementAt(i))); - featureList.push_back(std::make_pair("RunLength. (" + strRange + ") LongRunHighGreyLevelEmphasis Comb.", featureStd->ElementAt(i))); + featureList.push_back(std::make_pair("RunLength. (" + strRange + ") LongRunHighGreyLevelEmphasis Comb.", featureCombined->ElementAt(i))); break; case TextureFilterType::RunPercentage : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunPercentage Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") RunPercentage Std.", featureStd->ElementAt(i))); - featureList.push_back(std::make_pair("RunLength. (" + strRange + ") RunPercentage Comb.", featureStd->ElementAt(i))); + featureList.push_back(std::make_pair("RunLength. (" + strRange + ") RunPercentage Comb.", featureCombined->ElementAt(i)/newOffset->size())); break; case TextureFilterType::NumberOfRuns : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") NumberOfRuns Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") NumberOfRuns Std.", featureStd->ElementAt(i))); - featureList.push_back(std::make_pair("RunLength. (" + strRange + ") NumberOfRuns Comb.", featureStd->ElementAt(i))); + featureList.push_back(std::make_pair("RunLength. (" + strRange + ") NumberOfRuns Comb.", featureCombined->ElementAt(i))); break; case TextureFilterType::GreyLevelVariance : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") GreyLevelVariance Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") GreyLevelVariance Std.", featureStd->ElementAt(i))); - featureList.push_back(std::make_pair("RunLength. (" + strRange + ") GreyLevelVariance Comb.", featureStd->ElementAt(i))); + featureList.push_back(std::make_pair("RunLength. (" + strRange + ") GreyLevelVariance Comb.", featureCombined->ElementAt(i))); break; case TextureFilterType::RunLengthVariance : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunLengthVariance Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") RunLengthVariance Std.", featureStd->ElementAt(i))); - featureList.push_back(std::make_pair("RunLength. (" + strRange + ") RunLengthVariance Comb.", featureStd->ElementAt(i))); + featureList.push_back(std::make_pair("RunLength. (" + strRange + ") RunLengthVariance Comb.", featureCombined->ElementAt(i))); break; case TextureFilterType::RunEntropy : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunEntropy Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. (" + strRange + ") RunEntropy Std.", featureStd->ElementAt(i))); - featureList.push_back(std::make_pair("RunLength. (" + strRange + ") RunEntropy Comb.", featureStd->ElementAt(i))); + featureList.push_back(std::make_pair("RunLength. (" + strRange + ") RunEntropy Comb.", featureCombined->ElementAt(i))); break; default: break; } } } mitk::GIFGrayLevelRunLength::GIFGrayLevelRunLength(): m_Range(1.0), m_UseCtRange(false) { SetShortName("rl"); SetLongName("run-length"); } mitk::GIFGrayLevelRunLength::FeatureListType mitk::GIFGrayLevelRunLength::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; ParameterStruct params; params.m_UseCtRange=m_UseCtRange; params.m_Range = m_Range; params.m_Direction = GetDirection(); params.MinimumIntensity = GetMinimumIntensity(); params.MaximumIntensity = GetMaximumIntensity(); params.UseMinimumIntensity = GetUseMinimumIntensity(); params.UseMaximumIntensity = GetUseMaximumIntensity(); params.Bins = GetBins(); MITK_INFO << params.MinimumIntensity; MITK_INFO << params.MaximumIntensity; MITK_INFO << params.m_UseCtRange; MITK_INFO << params.m_Direction; MITK_INFO << params.Bins; - MITK_INFO< bins; + bins.push_back(GetBins()); if (parsedArgs.count(name + "::bins")) { bins = SplitDouble(parsedArgs[name + "::bins"].ToString(), ';'); } - else - { - bins.push_back(1); - } for (std::size_t i = 0; i < bins.size(); ++i) { MITK_INFO << "Start calculating Run-length with range " << bins[i] << "...."; this->SetBins(bins[i]); auto localResults = this->CalculateFeatures(feature, maskNoNAN); featureList.insert(featureList.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating Run-length with range " << bins[i] << "...."; } } } diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbouringGreyLevelDependenceFeatures.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbouringGreyLevelDependenceFeatures.cpp index 6d9a30f047..0b9679f634 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbouringGreyLevelDependenceFeatures.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbouringGreyLevelDependenceFeatures.cpp @@ -1,408 +1,412 @@ #include // MITK #include #include #include // ITK #include #include #include #include // STL #include static void MatrixFeaturesTo(mitk::NGLDMMatrixFeatures features, std::string prefix, mitk::GIFNeighbouringGreyLevelDependenceFeature::FeatureListType &featureList); mitk::NGLDMMatrixHolder::NGLDMMatrixHolder(double min, double max, int number, int depenence) : m_MinimumRange(min), m_MaximumRange(max), m_Stepsize(0), m_NumberOfDependences(depenence), m_NumberOfBins(number), m_NeighbourhoodSize(1), m_NumberOfNeighbourVoxels(0), m_NumberOfDependenceNeighbourVoxels(0), m_NumberOfNeighbourhoods(0), m_NumberOfCompleteNeighbourhoods(0) { m_Matrix.resize(number, depenence); m_Matrix.fill(0); m_Stepsize = (max - min) / (number); } int mitk::NGLDMMatrixHolder::IntensityToIndex(double intensity) { return std::floor((intensity - m_MinimumRange) / m_Stepsize); } double mitk::NGLDMMatrixHolder::IndexToMinIntensity(int index) { return m_MinimumRange + index * m_Stepsize; } double mitk::NGLDMMatrixHolder::IndexToMeanIntensity(int index) { return m_MinimumRange + (index+0.5) * m_Stepsize; } double mitk::NGLDMMatrixHolder::IndexToMaxIntensity(int index) { return m_MinimumRange + (index + 1) * m_Stepsize; } template void CalculateNGLDMMatrix(itk::Image* itkImage, itk::Image* mask, int alpha, int range, unsigned int direction, mitk::NGLDMMatrixHolder &holder) { typedef itk::Image ImageType; typedef itk::Image MaskImageType; typedef itk::NeighborhoodIterator ShapeIterType; typedef itk::NeighborhoodIterator ShapeMaskIterType; holder.m_NumberOfCompleteNeighbourhoods = 0; holder.m_NumberOfNeighbourhoods = 0; holder.m_NumberOfNeighbourVoxels = 0; holder.m_NumberOfDependenceNeighbourVoxels = 0; itk::Size radius; radius.Fill(range); if ((direction > 1) && (direction - 2 GetLargestPossibleRegion()); ShapeMaskIterType maskIter(radius, mask, mask->GetLargestPossibleRegion()); auto region = mask->GetLargestPossibleRegion(); auto center = imageIter.Size() / 2; auto iterSize = imageIter.Size(); holder.m_NeighbourhoodSize = iterSize-1; while (!maskIter.IsAtEnd()) { int sameValues = 0; bool completeNeighbourhood = true; int i = holder.IntensityToIndex(imageIter.GetCenterPixel()); if ((imageIter.GetCenterPixel() != imageIter.GetCenterPixel()) || (maskIter.GetCenterPixel() < 1)) { ++imageIter; ++maskIter; continue; } for (unsigned int position = 0; position < iterSize; ++position) { if (position == center) { continue; } if ( ! region.IsInside(maskIter.GetIndex(position))) { completeNeighbourhood = false; continue; } bool isInBounds; auto jIntensity = imageIter.GetPixel(position, isInBounds); auto jMask = maskIter.GetPixel(position, isInBounds); if (jMask < 1 || (jIntensity != jIntensity) || ( ! isInBounds)) { completeNeighbourhood = false; continue; } int j = holder.IntensityToIndex(jIntensity); holder.m_NumberOfNeighbourVoxels += 1; if (std::abs(i - j) <= alpha) { holder.m_NumberOfDependenceNeighbourVoxels += 1; ++sameValues; } } holder.m_Matrix(i, sameValues) += 1; holder.m_NumberOfNeighbourhoods += 1; if (completeNeighbourhood) { holder.m_NumberOfCompleteNeighbourhoods += 1; } ++imageIter; ++maskIter; } } void LocalCalculateFeatures( mitk::NGLDMMatrixHolder &holder, mitk::NGLDMMatrixFeatures & results ) { auto sijMatrix = holder.m_Matrix; auto piMatrix = holder.m_Matrix; auto pjMatrix = holder.m_Matrix; // double Ng = holder.m_NumberOfBins; // int NgSize = holder.m_NumberOfBins; double Ns = sijMatrix.sum(); piMatrix.rowwise().normalize(); pjMatrix.colwise().normalize(); for (int i = 0; i < holder.m_NumberOfBins; ++i) { double sj = 0; for (int j = 0; j < holder.m_NumberOfDependences; ++j) { - double iInt = holder.IndexToMeanIntensity(i); + double iInt = i + 1;// holder.IndexToMeanIntensity(i); double sij = sijMatrix(i, j); double k = j+1; double pij = sij / Ns; results.LowDependenceEmphasis += sij / k / k; results.HighDependenceEmphasis += sij * k*k; if (iInt != 0) { results.LowGreyLevelCountEmphasis += sij / iInt / iInt; } results.HighGreyLevelCountEmphasis += sij * iInt*iInt; if (iInt != 0) { results.LowDependenceLowGreyLevelEmphasis += sij / k / k / iInt / iInt; } results.LowDependenceHighGreyLevelEmphasis += sij * iInt*iInt / k / k; if (iInt != 0) { results.HighDependenceLowGreyLevelEmphasis += sij *k * k / iInt / iInt; } results.HighDependenceHighGreyLevelEmphasis += sij * k*k*iInt*iInt; results.MeanGreyLevelCount += iInt * pij; results.MeanDependenceCount += k * pij; if (pij > 0) { results.DependenceCountEntropy -= pij * std::log(pij) / std::log(2); } results.DependenceCountEnergy += pij*pij; sj += sij; } results.GreyLevelNonUniformity += sj*sj; results.GreyLevelNonUniformityNormalised += sj*sj; } for (int j = 0; j < holder.m_NumberOfDependences; ++j) { double si = 0; for (int i = 0; i < holder.m_NumberOfBins; ++i) { double sij = sijMatrix(i, j); si += sij; } results.DependenceCountNonUniformity += si*si; results.DependenceCountNonUniformityNormalised += si*si; } for (int i = 0; i < holder.m_NumberOfBins; ++i) { for (int j = 0; j < holder.m_NumberOfDependences; ++j) { - double iInt = holder.IndexToMeanIntensity(i); + double iInt = i + 1;// holder.IndexToMeanIntensity(i); double sij = sijMatrix(i, j); double k = j + 1; double pij = sij / Ns; results.GreyLevelVariance += (iInt - results.MeanGreyLevelCount)* (iInt - results.MeanGreyLevelCount) * pij; results.DependenceCountVariance += (k - results.MeanDependenceCount)* (k - results.MeanDependenceCount) * pij; } } results.LowDependenceEmphasis /= Ns; results.HighDependenceEmphasis /= Ns; results.LowGreyLevelCountEmphasis /= Ns; results.HighGreyLevelCountEmphasis /= Ns; results.LowDependenceLowGreyLevelEmphasis /= Ns; results.LowDependenceHighGreyLevelEmphasis /= Ns; results.HighDependenceLowGreyLevelEmphasis /= Ns; results.HighDependenceHighGreyLevelEmphasis /= Ns; results.GreyLevelNonUniformity /= Ns; results.GreyLevelNonUniformityNormalised /= (Ns*Ns); results.DependenceCountNonUniformity /= Ns; results.DependenceCountNonUniformityNormalised /= (Ns*Ns); results.DependenceCountPercentage = 1; results.ExpectedNeighbourhoodSize = holder.m_NeighbourhoodSize; results.AverageNeighbourhoodSize = holder.m_NumberOfNeighbourVoxels / (1.0 * holder.m_NumberOfNeighbourhoods); results.AverageIncompleteNeighbourhoodSize = (holder.m_NumberOfNeighbourVoxels - holder.m_NumberOfCompleteNeighbourhoods* holder.m_NeighbourhoodSize) / (1.0 * (holder.m_NumberOfNeighbourhoods - holder.m_NumberOfCompleteNeighbourhoods)); results.PercentageOfCompleteNeighbourhoods = (1.0*holder.m_NumberOfCompleteNeighbourhoods) / (1.0 * holder.m_NumberOfNeighbourhoods); results.PercentageOfDependenceNeighbours = holder.m_NumberOfDependenceNeighbourVoxels / (1.0 * holder.m_NumberOfNeighbourVoxels); } template void CalculateCoocurenceFeatures(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFNeighbouringGreyLevelDependenceFeature::FeatureListType & featureList, mitk::GIFNeighbouringGreyLevelDependenceFeature::GIFNeighbouringGreyLevelDependenceFeatureConfiguration config) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::MinimumMaximumImageCalculator MinMaxComputerType; /////////////////////////////////////////////////////////////////////////////////////////////// typename MinMaxComputerType::Pointer minMaxComputer = MinMaxComputerType::New(); minMaxComputer->SetImage(itkImage); minMaxComputer->Compute(); double rangeMin = -0.5 + minMaxComputer->GetMinimum(); double rangeMax = 0.5 + minMaxComputer->GetMaximum(); if (config.UseMinimumIntensity) rangeMin = config.MinimumIntensity; if (config.UseMaximumIntensity) rangeMax = config.MaximumIntensity; // Define Range int numberOfBins = config.Bins; typename MaskType::Pointer maskImage = MaskType::New(); mitk::CastToItkImage(mask, maskImage); std::vector resultVector; + int numberofDependency = 37; + if (VImageDimension == 2) + numberofDependency = 37; + mitk::NGLDMMatrixHolder holderOverall(rangeMin, rangeMax, numberOfBins,37); mitk::NGLDMMatrixFeatures overallFeature; CalculateNGLDMMatrix(itkImage, maskImage, config.alpha, config.range, config.direction, holderOverall); LocalCalculateFeatures(holderOverall, overallFeature); std::ostringstream ss; ss << config.range; std::string strRange = ss.str(); MatrixFeaturesTo(overallFeature, "NGLDM (" + strRange + ") overall", featureList); } static void MatrixFeaturesTo(mitk::NGLDMMatrixFeatures features, std::string prefix, mitk::GIFNeighbouringGreyLevelDependenceFeature::FeatureListType &featureList) { featureList.push_back(std::make_pair(prefix + " Low Dependence Emphasis", features.LowDependenceEmphasis)); featureList.push_back(std::make_pair(prefix + " High Dependence Emphasis", features.HighDependenceEmphasis)); featureList.push_back(std::make_pair(prefix + " Low Grey Level Count Emphasis", features.LowGreyLevelCountEmphasis)); featureList.push_back(std::make_pair(prefix + " High Grey Level Count Emphasis", features.HighGreyLevelCountEmphasis)); featureList.push_back(std::make_pair(prefix + " Low Dependence Low Grey Level Emphasis", features.LowDependenceLowGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + " Low Dependence High Grey Level Emphasis", features.LowDependenceHighGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + " High Dependence Low Grey Level Emphasis", features.HighDependenceLowGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + " High Dependence High Grey Level Emphasis", features.HighDependenceHighGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + " Grey Level Non-Uniformity", features.GreyLevelNonUniformity)); featureList.push_back(std::make_pair(prefix + " Grey Level Non-Uniformity Normalised", features.GreyLevelNonUniformityNormalised)); featureList.push_back(std::make_pair(prefix + " Dependence Count Non-Uniformity", features.DependenceCountNonUniformity)); featureList.push_back(std::make_pair(prefix + " Dependence Count Non-Uniformtiy Normalised", features.DependenceCountNonUniformityNormalised)); featureList.push_back(std::make_pair(prefix + " Dependence Count Percentage", features.DependenceCountPercentage)); featureList.push_back(std::make_pair(prefix + " Grey Level Mean", features.MeanGreyLevelCount)); featureList.push_back(std::make_pair(prefix + " Grey Level Variance", features.GreyLevelVariance)); featureList.push_back(std::make_pair(prefix + " Dependence Count Mean", features.MeanDependenceCount)); featureList.push_back(std::make_pair(prefix + " Dependence Count Variance", features.DependenceCountVariance)); featureList.push_back(std::make_pair(prefix + " Dependence Count Entropy", features.DependenceCountEntropy)); featureList.push_back(std::make_pair(prefix + " Dependence Count Energy", features.DependenceCountEnergy)); featureList.push_back(std::make_pair(prefix + " Expected Neighbourhood Size", features.ExpectedNeighbourhoodSize)); featureList.push_back(std::make_pair(prefix + " Average Neighbourhood Size", features.AverageNeighbourhoodSize)); featureList.push_back(std::make_pair(prefix + " Average Incomplete Neighbourhood Size", features.AverageIncompleteNeighbourhoodSize)); featureList.push_back(std::make_pair(prefix + " Percentage of complete Neighbourhoods", features.PercentageOfCompleteNeighbourhoods)); featureList.push_back(std::make_pair(prefix + " Percentage of Dependence Neighbour Voxels", features.PercentageOfDependenceNeighbours)); } mitk::GIFNeighbouringGreyLevelDependenceFeature::GIFNeighbouringGreyLevelDependenceFeature() : m_Range(1.0), m_Bins(6) { SetShortName("ngldm"); SetLongName("ngldm"); } mitk::GIFNeighbouringGreyLevelDependenceFeature::FeatureListType mitk::GIFNeighbouringGreyLevelDependenceFeature::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; GIFNeighbouringGreyLevelDependenceFeatureConfiguration config; config.direction = GetDirection(); config.range = m_Range; config.alpha = 0; config.MinimumIntensity = GetMinimumIntensity(); config.MaximumIntensity = GetMaximumIntensity(); config.UseMinimumIntensity = GetUseMinimumIntensity(); config.UseMaximumIntensity = GetUseMaximumIntensity(); config.Bins = GetBins(); AccessByItk_3(image, CalculateCoocurenceFeatures, mask, featureList,config); return featureList; } mitk::GIFNeighbouringGreyLevelDependenceFeature::FeatureNameListType mitk::GIFNeighbouringGreyLevelDependenceFeature::GetFeatureNames() { FeatureNameListType featureList; return featureList; } void mitk::GIFNeighbouringGreyLevelDependenceFeature::AddArguments(mitkCommandLineParser &parser) { std::string name = GetOptionPrefix(); parser.addArgument(GetLongName(), name, mitkCommandLineParser::String, "Calculate Neighbouring grey level dependence based features", "Calculate Neighbouring grey level dependence based features", us::Any()); parser.addArgument(name + "::range", name + "::range", mitkCommandLineParser::String, "NGLD Range", "Define the range that is used (Semicolon-separated)", us::Any()); parser.addArgument(name + "::bins", name + "::bins", mitkCommandLineParser::String, "NGLD Number of Bins", "Define the number of bins that is used ", us::Any()); } void mitk::GIFNeighbouringGreyLevelDependenceFeature::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &, const Image::Pointer &maskNoNAN, FeatureListType &featureList) { auto parsedArgs = GetParameter(); std::string name = GetOptionPrefix(); if (parsedArgs.count(GetLongName())) { std::vector ranges; if (parsedArgs.count(name + "::range")) { ranges = SplitDouble(parsedArgs[name + "::range"].ToString(), ';'); } else { ranges.push_back(1); } if (parsedArgs.count(name + "::bins")) { auto bins = SplitDouble(parsedArgs[name + "::bins"].ToString(), ';')[0]; this->SetBins(bins); } for (std::size_t i = 0; i < ranges.size(); ++i) { MITK_INFO << "Start calculating NGLDM with range " << ranges[i] << "...."; this->SetRange(ranges[i]); auto localResults = this->CalculateFeatures(feature, maskNoNAN); featureList.insert(featureList.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating NGLDM with range " << ranges[i] << "...."; } } } diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricStatistics.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricStatistics.cpp index 80b47d5fef..ca5faec970 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricStatistics.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricStatistics.cpp @@ -1,432 +1,432 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include // MITK #include #include #include #include #include // ITK #include #include // VTK #include #include #include // STL #include #include // OpenCV #include template void CalculateVolumeStatistic(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFVolumetricStatistics::FeatureListType & featureList) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::LabelStatisticsImageFilter FilterType; typename MaskType::Pointer maskImage = MaskType::New(); mitk::CastToItkImage(mask, maskImage); typename FilterType::Pointer labelStatisticsImageFilter = FilterType::New(); labelStatisticsImageFilter->SetInput( itkImage ); labelStatisticsImageFilter->SetLabelInput(maskImage); labelStatisticsImageFilter->Update(); double volume = labelStatisticsImageFilter->GetCount(1); double voxelVolume = 1; for (int i = 0; i < (int)(VImageDimension); ++i) { volume *= itkImage->GetSpacing()[i]; voxelVolume *= itkImage->GetSpacing()[i]; } featureList.push_back(std::make_pair("Volumetric Features Volume (pixel based)", volume)); featureList.push_back(std::make_pair("Volumetric Features Voxel Volume", voxelVolume)); } template void CalculateLargestDiameter(itk::Image* mask, mitk::Image::Pointer valueImage, mitk::GIFVolumetricStatistics::FeatureListType & featureList) { typedef itk::Image ValueImageType; typename ValueImageType::Pointer itkValueImage = ValueImageType::New(); mitk::CastToItkImage(valueImage, itkValueImage); typedef itk::Image ImageType; typedef typename ImageType::PointType PointType; typename ImageType::SizeType radius; for (int i=0; i < (int)VImageDimension; ++i) radius[i] = 1; itk::NeighborhoodIterator iterator(radius, mask, mask->GetRequestedRegion()); itk::NeighborhoodIterator valueIter(radius, itkValueImage, itkValueImage->GetRequestedRegion()); std::vector borderPoints; // // Calculate surface in different directions // double surface = 0; std::vector directionSurface; for (int i = 0; i < (int)(iterator.Size()); ++i) { auto offset = iterator.GetOffset(i); double deltaS = 1; int nonZeros = 0; for (unsigned int j = 0; j < VImageDimension; ++j) { if (offset[j] != 0 && nonZeros == 0) { for (unsigned int k = 0; k < VImageDimension; ++k) { if (k != j) deltaS *= mask->GetSpacing()[k]; } nonZeros++; } else if (offset[j] != 0) { deltaS = 0; } } if (nonZeros < 1) deltaS = 0; directionSurface.push_back(deltaS); } // // Prepare calulation of Centre of mass shift // PointType normalCenter(0); PointType normalCenterUncorrected(0); PointType weightedCenter(0); PointType currentPoint; int numberOfPoints = 0; int numberOfPointsUncorrected = 0; double sumOfPoints = 0; while(!iterator.IsAtEnd()) { if (iterator.GetCenterPixel() == 0) { ++iterator; ++valueIter; continue; } mask->TransformIndexToPhysicalPoint(iterator.GetIndex(), currentPoint); normalCenterUncorrected += currentPoint.GetVectorFromOrigin(); ++numberOfPointsUncorrected; double intensityValue = valueIter.GetCenterPixel(); if (intensityValue == intensityValue) { normalCenter += currentPoint.GetVectorFromOrigin(); weightedCenter += currentPoint.GetVectorFromOrigin() * intensityValue; sumOfPoints += intensityValue; ++numberOfPoints; } bool border = false; for (int i = 0; i < (int)(iterator.Size()); ++i) { if (iterator.GetPixel(i) == 0 || ( ! iterator.IndexInBounds(i))) { border = true; surface += directionSurface[i]; //break; } } if (border) { auto centerIndex = iterator.GetIndex(); PointType centerPoint; mask->TransformIndexToPhysicalPoint(centerIndex, centerPoint ); borderPoints.push_back(centerPoint); } ++iterator; ++valueIter; } auto normalCenterVector = normalCenter.GetVectorFromOrigin() / numberOfPoints; auto normalCenterVectorUncorrected = normalCenter.GetVectorFromOrigin() / numberOfPointsUncorrected; auto weightedCenterVector = weightedCenter.GetVectorFromOrigin() / sumOfPoints; auto differenceOfCentersUncorrected = (normalCenterVectorUncorrected - weightedCenterVector).GetNorm(); auto differenceOfCenters = (normalCenterVector - weightedCenterVector).GetNorm(); double longestDiameter = 0; unsigned long numberOfBorderPoints = borderPoints.size(); for (int i = 0; i < (int)numberOfBorderPoints; ++i) { auto point = borderPoints[i]; for (int j = i; j < (int)numberOfBorderPoints; ++j) { double newDiameter=point.EuclideanDistanceTo(borderPoints[j]); if (newDiameter > longestDiameter) longestDiameter = newDiameter; } } featureList.push_back(std::make_pair("Volumetric Features Maximum 3D diameter", longestDiameter)); featureList.push_back(std::make_pair("Volumetric Features Surface (Voxel based)", surface)); featureList.push_back(std::make_pair("Volumetric Features Centre of mass shift", differenceOfCenters)); featureList.push_back(std::make_pair("Volumetric Features Centre of mass shift (Uncorrected)", differenceOfCentersUncorrected)); } mitk::GIFVolumetricStatistics::GIFVolumetricStatistics() { SetLongName("volume"); SetShortName("vol"); } mitk::GIFVolumetricStatistics::FeatureListType mitk::GIFVolumetricStatistics::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; if (image->GetDimension() < 3) { return featureList; } AccessByItk_2(image, CalculateVolumeStatistic, mask, featureList); AccessByItk_2(mask, CalculateLargestDiameter, image, featureList); vtkSmartPointer mesher = vtkSmartPointer::New(); vtkSmartPointer stats = vtkSmartPointer::New(); mesher->SetInputData(mask->GetVtkImageData()); mesher->SetValue(0, 0.5); stats->SetInputConnection(mesher->GetOutputPort()); stats->Update(); double pi = vnl_math::pi; double meshVolume = stats->GetVolume(); double meshSurf = stats->GetSurfaceArea(); double pixelVolume = featureList[0].second; double pixelSurface = featureList[3].second; MITK_INFO << "Surf: " << pixelSurface << " Vol " << pixelVolume; double compactness1 = pixelVolume / (std::sqrt(pi) * std::pow(meshSurf, 2.0 / 3.0)); double compactness1Pixel = pixelVolume / (std::sqrt(pi) * std::pow(pixelSurface, 2.0 / 3.0)); //This is the definition used by Aertz. However, due to 2/3 this feature is not demensionless. Use compactness3 instead. double compactness2 = 36 * pi*pixelVolume*pixelVolume / meshSurf / meshSurf / meshSurf; double compactness2Pixel = 36 * pi*pixelVolume*pixelVolume / pixelSurface / pixelSurface / pixelSurface; double compactness3 = pixelVolume / (std::sqrt(pi) * std::pow(meshSurf, 3.0 / 2.0)); double compactness3Pixel = pixelVolume / (std::sqrt(pi) * std::pow(pixelSurface, 3.0 / 2.0)); double sphericity = std::pow(pi, 1 / 3.0) *std::pow(6 * pixelVolume, 2.0 / 3.0) / meshSurf; double sphericityPixel = std::pow(pi, 1 / 3.0) *std::pow(6 * pixelVolume, 2.0 / 3.0) / pixelSurface; double surfaceToVolume = meshSurf / pixelVolume; double surfaceToVolumePixel = pixelSurface / pixelVolume; double sphericalDisproportion = meshSurf / 4 / pi / std::pow(3.0 / 4.0 / pi * pixelVolume, 2.0 / 3.0); double sphericalDisproportionPixel = pixelSurface / 4 / pi / std::pow(3.0 / 4.0 / pi * pixelVolume, 2.0 / 3.0); double asphericity = std::pow(1.0/compactness2, (1.0 / 3.0)) - 1; double asphericityPixel = std::pow(1.0/compactness2Pixel, (1.0 / 3.0)) - 1; //Calculate center of mass shift int xx = mask->GetDimensions()[0]; int yy = mask->GetDimensions()[1]; int zz = mask->GetDimensions()[2]; double xd = mask->GetGeometry()->GetSpacing()[0]; double yd = mask->GetGeometry()->GetSpacing()[1]; double zd = mask->GetGeometry()->GetSpacing()[2]; std::vector pointsForPCA; std::vector pointsForPCAUncorrected; for (int x = 0; x < xx; x++) { for (int y = 0; y < yy; y++) { for (int z = 0; z < zz; z++) { itk::Image::IndexType index; index[0] = x; index[1] = y; index[2] = z; mitk::ScalarType pxImage; mitk::ScalarType pxMask; mitkPixelTypeMultiplex5( mitk::FastSinglePixelAccess, image->GetChannelDescriptor().GetPixelType(), image, image->GetVolumeData(), index, pxImage, 0); mitkPixelTypeMultiplex5( mitk::FastSinglePixelAccess, mask->GetChannelDescriptor().GetPixelType(), mask, mask->GetVolumeData(), index, pxMask, 0); //Check if voxel is contained in segmentation if (pxMask > 0) { cv::Point3d tmp; tmp.x = x * xd; tmp.y = y * yd; tmp.z = z * zd; pointsForPCAUncorrected.push_back(tmp); if (pxImage == pxImage) { pointsForPCA.push_back(tmp); } } } } } //Calculate PCA Features int sz = pointsForPCA.size(); cv::Mat data_pts = cv::Mat(sz, 3, CV_64FC1); for (int i = 0; i < data_pts.rows; ++i) { data_pts.at(i, 0) = pointsForPCA[i].x; data_pts.at(i, 1) = pointsForPCA[i].y; data_pts.at(i, 2) = pointsForPCA[i].z; } //Calculate PCA Features int szUC = pointsForPCAUncorrected.size(); cv::Mat data_ptsUC = cv::Mat(szUC, 3, CV_64FC1); for (int i = 0; i < data_ptsUC.rows; ++i) { data_ptsUC.at(i, 0) = pointsForPCAUncorrected[i].x; data_ptsUC.at(i, 1) = pointsForPCAUncorrected[i].y; data_ptsUC.at(i, 2) = pointsForPCAUncorrected[i].z; } //Perform PCA analysis cv::PCA pca_analysis(data_pts, cv::Mat(), CV_PCA_DATA_AS_ROW); cv::PCA pca_analysisUC(data_ptsUC, cv::Mat(), CV_PCA_DATA_AS_ROW); //Store the eigenvalues std::vector eigen_val(3); std::vector eigen_valUC(3); for (int i = 0; i < 3; ++i) { eigen_val[i] = pca_analysis.eigenvalues.at(0, i); eigen_valUC[i] = pca_analysisUC.eigenvalues.at(0, i); } std::sort(eigen_val.begin(), eigen_val.end()); std::sort(eigen_valUC.begin(), eigen_valUC.end()); - double major = eigen_val[2]; - double minor = eigen_val[1]; - double least = eigen_val[0]; - double elongation = major == 0 ? 0 : minor/major; - double flatness = major == 0 ? 0 : least / major; - double majorUC = eigen_valUC[2]; - double minorUC = eigen_valUC[1]; - double leastUC = eigen_valUC[0]; - double elongationUC = majorUC == 0 ? 0 : minorUC / majorUC; - double flatnessUC = majorUC == 0 ? 0 : leastUC / majorUC; + double major = 4*sqrt(eigen_val[2]); + double minor = 4*sqrt(eigen_val[1]); + double least = 4*sqrt(eigen_val[0]); + double elongation = major == 0 ? 0 : sqrt(minor/major); + double flatness = major == 0 ? 0 : sqrt(least / major); + double majorUC = 4*sqrt(eigen_valUC[2]); + double minorUC = 4*sqrt(eigen_valUC[1]); + double leastUC = 4*sqrt(eigen_valUC[0]); + double elongationUC = majorUC == 0 ? 0 : sqrt(minorUC / majorUC); + double flatnessUC = majorUC == 0 ? 0 : sqrt(leastUC / majorUC); featureList.push_back(std::make_pair("Volumetric Features Volume (mesh based)",meshVolume)); featureList.push_back(std::make_pair("Volumetric Features Surface area",meshSurf)); featureList.push_back(std::make_pair("Volumetric Features Surface to volume ratio",surfaceToVolume)); featureList.push_back(std::make_pair("Volumetric Features Sphericity",sphericity)); featureList.push_back(std::make_pair("Volumetric Features Asphericity",asphericity)); featureList.push_back(std::make_pair("Volumetric Features Compactness 1",compactness1)); featureList.push_back(std::make_pair("Volumetric Features Compactness 2",compactness2)); featureList.push_back(std::make_pair("Volumetric Features Compactness 3",compactness3)); featureList.push_back(std::make_pair("Volumetric Features Spherical disproportion", sphericalDisproportion)); featureList.push_back(std::make_pair("Volumetric Features Surface to volume ratio (Voxel based)", surfaceToVolumePixel)); featureList.push_back(std::make_pair("Volumetric Features Sphericity (Voxel based)", sphericityPixel)); featureList.push_back(std::make_pair("Volumetric Features Asphericity (Voxel based)", asphericityPixel)); featureList.push_back(std::make_pair("Volumetric Features Compactness 1 (Voxel based)", compactness1Pixel)); featureList.push_back(std::make_pair("Volumetric Features Compactness 2 (Voxel based)", compactness2Pixel)); featureList.push_back(std::make_pair("Volumetric Features Compactness 3 (Voxel based)", compactness3Pixel)); featureList.push_back(std::make_pair("Volumetric Features Spherical disproportion (Voxel based)", sphericalDisproportionPixel)); featureList.push_back(std::make_pair("Volumetric Features PCA Major Axis",major)); featureList.push_back(std::make_pair("Volumetric Features PCA Minor Axis",minor)); featureList.push_back(std::make_pair("Volumetric Features PCA Least Axis",least)); featureList.push_back(std::make_pair("Volumetric Features PCA Elongation",elongation)); featureList.push_back(std::make_pair("Volumetric Features PCA Flatness",flatness)); featureList.push_back(std::make_pair("Volumetric Features PCA Major Axis (Uncorrected)", majorUC)); featureList.push_back(std::make_pair("Volumetric Features PCA Minor Axis (Uncorrected)", minorUC)); featureList.push_back(std::make_pair("Volumetric Features PCA Least Axis (Uncorrected)", leastUC)); featureList.push_back(std::make_pair("Volumetric Features PCA Elongation (Uncorrected)", elongationUC)); featureList.push_back(std::make_pair("Volumetric Features PCA Flatness (Uncorrected)", flatnessUC)); return featureList; } mitk::GIFVolumetricStatistics::FeatureNameListType mitk::GIFVolumetricStatistics::GetFeatureNames() { FeatureNameListType featureList; featureList.push_back("Volumetric Features Compactness 1"); featureList.push_back("Volumetric Features Compactness 2"); featureList.push_back("Volumetric Features Compactness 3"); featureList.push_back("Volumetric Features Sphericity"); featureList.push_back("Volumetric Features Asphericity"); featureList.push_back("Volumetric Features Surface to volume ratio"); featureList.push_back("Volumetric Features Surface area"); featureList.push_back("Volumetric Features Volume (mesh based)"); featureList.push_back("Volumetric Features Volume (pixel based)"); featureList.push_back("Volumetric Features Spherical disproportion"); featureList.push_back("Volumetric Features Maximum 3D diameter"); return featureList; } void mitk::GIFVolumetricStatistics::AddArguments(mitkCommandLineParser &parser) { std::string name = GetOptionPrefix(); parser.addArgument(GetLongName(), name, mitkCommandLineParser::String, "Use Volume-Statistic", "calculates volume based features", us::Any()); } void mitk::GIFVolumetricStatistics::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &, FeatureListType &featureList) { auto parsedArgs = GetParameter(); if (parsedArgs.count(GetLongName())) { MITK_INFO << "Start calculating volumetric features ...."; auto localResults = this->CalculateFeatures(feature, mask); featureList.insert(featureList.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating volumetric features...."; } }