X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=tools%2FclitkImageStatisticsGenericFilter.txx;h=d1b0d4c4abf845ce98e50054e8d259093523e37a;hb=c2e0628b1d9f0940ac192ff3683c3ede5d01ceb3;hp=46abbc2538cc959881b53f55d15ff722a1fffb75;hpb=a26cd8a19e1b9ad8344ab501436045f171a73713;p=clitk.git diff --git a/tools/clitkImageStatisticsGenericFilter.txx b/tools/clitkImageStatisticsGenericFilter.txx old mode 100755 new mode 100644 index 46abbc2..d1b0d4c --- a/tools/clitkImageStatisticsGenericFilter.txx +++ b/tools/clitkImageStatisticsGenericFilter.txx @@ -3,7 +3,7 @@ Authors belong to: - University of LYON http://www.universite-lyon.fr/ - - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr + - Léon Bérard cancer center http://www.centreleonberard.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even @@ -14,10 +14,13 @@ - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html -======================================================================-====*/ +===========================================================================**/ #ifndef clitkImageStatisticsGenericFilter_txx #define clitkImageStatisticsGenericFilter_txx +#include "itkNthElementImageAdaptor.h" +#include "itkJoinSeriesImageFilter.h" + /* ================================================= * @file clitkImageStatisticsGenericFilter.txx * @author @@ -26,6 +29,7 @@ * @brief * ===================================================*/ +#include "clitkImageStatisticsGenericFilter.h" namespace clitk @@ -34,15 +38,15 @@ namespace clitk //------------------------------------------------------------------- // Update with the number of dimensions //------------------------------------------------------------------- - template + template void ImageStatisticsGenericFilter::UpdateWithDim(std::string PixelType) { - if (m_Verbose) std::cout << "Image was detected to be "<(); + UpdateWithDimAndPixelType(); } // else if(PixelType == "unsigned_short"){ // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl; @@ -51,7 +55,7 @@ namespace clitk else if (PixelType == "unsigned_char"){ if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl; - UpdateWithDimAndPixelType(); + UpdateWithDimAndPixelType(); } // else if (PixelType == "char"){ @@ -60,7 +64,7 @@ namespace clitk // } else { if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl; - UpdateWithDimAndPixelType(); + UpdateWithDimAndPixelType(); } } @@ -68,15 +72,15 @@ namespace clitk //------------------------------------------------------------------- // Update with the number of dimensions and the pixeltype //------------------------------------------------------------------- - template + template void ImageStatisticsGenericFilter::UpdateWithDimAndPixelType() { // ImageTypes - typedef itk::Image InputImageType; - typedef itk::Image LabelImageType; - typedef itk::Image OutputImageType; + typedef unsigned char LabelPixelType; + typedef itk::Image, Dimension> InputImageType; + typedef itk::Image LabelImageType; // Read the input typedef itk::ImageFileReader InputReaderType; @@ -84,31 +88,64 @@ namespace clitk reader->SetFileName( m_InputFileName); reader->Update(); typename InputImageType::Pointer input= reader->GetOutput(); + + typedef itk::NthElementImageAdaptor InputImageAdaptorType; + typedef itk::Image OutputImageType; + typename InputImageAdaptorType::Pointer input_adaptor = InputImageAdaptorType::New(); + input_adaptor->SetImage(input); + // Filter - typedef itk::LabelStatisticsImageFilter StatisticsImageFilterType; + typedef itk::LabelStatisticsImageFilter StatisticsImageFilterType; typename StatisticsImageFilterType::Pointer statisticsFilter=StatisticsImageFilterType::New(); - statisticsFilter->SetInput(input); + statisticsFilter->SetInput(input_adaptor); // Label image typename LabelImageType::Pointer labelImage; - if (m_ArgsInfo.mask_given) - { - typedef itk::ImageFileReader LabelImageReaderType; - typename LabelImageReaderType::Pointer labelImageReader=LabelImageReaderType::New(); - labelImageReader->SetFileName(m_ArgsInfo.mask_arg); - labelImageReader->Update(); - labelImage= labelImageReader->GetOutput(); + if (m_ArgsInfo.mask_given) { + int maskDimension, maskComponents; + std::string maskPixelType; + ReadImageDimensionAndPixelType(m_ArgsInfo.mask_arg, maskDimension, maskPixelType, maskComponents); + if (maskDimension == Dimension - 1) { + // Due to a limitation of filter itk::LabelStatisticsImageFilter, InputImageType and LabelImageType + // must have the same image dimension. However, we want to support label images with Dl = Di - 1, + // so we need to replicate the label image as many times as the size along dimension Di. + if (m_Verbose) + std::cout << "Replicating label image to match input image's dimension... " << std::endl; + + typedef itk::Image ReducedLabelImageType; + typedef itk::ImageFileReader LabelImageReaderType; + typedef itk::JoinSeriesImageFilter JoinImageFilterType; + + typename LabelImageReaderType::Pointer labelImageReader=LabelImageReaderType::New(); + labelImageReader->SetFileName(m_ArgsInfo.mask_arg); + labelImageReader->Update(); + + typename JoinImageFilterType::Pointer joinFilter = JoinImageFilterType::New(); + typename InputImageType::SizeType size = input->GetLargestPossibleRegion().GetSize(); + for (unsigned int i = 0; i < size[Dimension - 1]; i++) + joinFilter->PushBackInput(labelImageReader->GetOutput()); + + joinFilter->Update(); + labelImage = joinFilter->GetOutput(); } - else - { - labelImage=LabelImageType::New(); - labelImage->SetRegions(input->GetLargestPossibleRegion()); - labelImage->SetOrigin(input->GetOrigin()); - labelImage->SetSpacing(input->GetSpacing()); - labelImage->Allocate(); - labelImage->FillBuffer(m_ArgsInfo.label_arg[0]); + else { + typedef itk::ImageFileReader LabelImageReaderType; + typename LabelImageReaderType::Pointer labelImageReader=LabelImageReaderType::New(); + labelImageReader->SetFileName(m_ArgsInfo.mask_arg); + labelImageReader->Update(); + labelImage= labelImageReader->GetOutput(); } + + } + else { + labelImage=LabelImageType::New(); + labelImage->SetRegions(input->GetLargestPossibleRegion()); + labelImage->SetOrigin(input->GetOrigin()); + labelImage->SetSpacing(input->GetSpacing()); + labelImage->Allocate(); + labelImage->FillBuffer(m_ArgsInfo.label_arg[0]); + } statisticsFilter->SetLabelInput(labelImage); // For each Label @@ -119,77 +156,91 @@ namespace clitk else numberOfLabels=1; - for (unsigned int k=0; k< numberOfLabels; k++) - { - label=m_ArgsInfo.label_arg[k]; - - std::cout<SetUseHistograms(true); - statisticsFilter->SetHistogramParameters(m_ArgsInfo.bins_arg, m_ArgsInfo.lower_arg, m_ArgsInfo.upper_arg); - } - statisticsFilter->Update(); - - - // Output - if (m_Verbose) std::cout<<"N° of pixels: "; - std::cout<GetCount(label)<GetMean(label)<GetSigma(label)<GetVariance(label)<GetMinimum(label)<GetMaximum(label)<GetSum(label)<GetBoundingBox(label).size(); i++) - std::cout<GetBoundingBox(label)[i]<<" "; - std::cout<GetMedian(label)<GetHistogram(label); - - // Screen - if (m_Verbose) std::cout<<"Histogram: "<GetBinMin(0,i)<<"\t"<GetMeasurement(i,0)<<"\t"<GetBinMax(0,i)<<"\t"<GetFrequency(i)<GetBinMin(0,i)<<"\t"<GetMeasurement(i,0)<<"\t"<GetBinMax(0,i)<<"\t"<GetFrequency(i)<SelectNthElement(c); + input_adaptor->Update(); + + for (unsigned int k=0; k< numberOfLabels; k++) { + label=m_ArgsInfo.label_arg[k]; + + std::cout<SetUseHistograms(true); + statisticsFilter->SetHistogramParameters(m_ArgsInfo.bins_arg, m_ArgsInfo.lower_arg, m_ArgsInfo.upper_arg); + } + statisticsFilter->Update(); + + // Output + if (m_Verbose) std::cout<<"N° of pixels: "; + std::cout<GetCount(label)<GetMean(label)<GetSigma(label)<GetVariance(label)<GetMinimum(label)<GetMaximum(label)<GetSum(label)<GetBoundingBox(label).size(); i++) + std::cout<GetBoundingBox(label)[i]<<" "; + std::cout<GetMedian(label)<GetHistogram(label); + + // Screen + if (m_Verbose) std::cout<<"Histogram: "<GetBinMin(0,i)<<"\t"<GetMeasurement(i,0)<<"\t"<GetBinMax(0,i)<<"\t"<GetFrequency(i)<GetBinMin(0,i)<<"\t"<GetMeasurement(i,0)<<"\t"<GetBinMax(0,i)<<"\t"<GetFrequency(i)<