1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://www.centreleonberard.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
9 This software is distributed WITHOUT ANY WARRANTY; without even
10 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 PURPOSE. See the copyright notices for more information.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================**/
19 #ifndef clitkMaskOfIntegratedIntensityGenericFilter_txx
20 #define clitkMaskOfIntegratedIntensityGenericFilter_txx
23 #include "itkIntensityWindowingImageFilter.h"
24 #include "itkLabelStatisticsImageFilter.h"
25 #include "itkMaskImageFilter.h"
26 #include "itkMaskNegatedImageFilter.h"
27 #include <clitkCommon.h>
33 //--------------------------------------------------------------------
34 template<class args_info_type>
35 MaskOfIntegratedIntensityGenericFilter<args_info_type>::MaskOfIntegratedIntensityGenericFilter():
36 ImageToImageGenericFilter<Self>("MaskOfIntegratedIntensity")
38 InitializeImageType<2>();
39 InitializeImageType<3>();
41 //--------------------------------------------------------------------
44 //--------------------------------------------------------------------
45 template<class args_info_type>
46 template<unsigned int Dim>
47 void MaskOfIntegratedIntensityGenericFilter<args_info_type>::InitializeImageType()
49 ADD_DEFAULT_IMAGE_TYPES(Dim);
51 //--------------------------------------------------------------------
54 //--------------------------------------------------------------------
55 template<class args_info_type>
56 void MaskOfIntegratedIntensityGenericFilter<args_info_type>::SetArgsInfo(const args_info_type & a)
59 this->SetIOVerbose(mArgsInfo.verbose_flag);
60 if (mArgsInfo.imagetypes_flag) this->PrintAvailableImageTypes();
62 this->SetInputFilename(mArgsInfo.input_arg);
63 this->SetOutputFilename(mArgsInfo.output_arg);
64 this->SetPercentage(mArgsInfo.percentage_arg);
67 //--------------------------------------------------------------------
70 //--------------------------------------------------------------------
71 // https://stackoverflow.com/questions/1577475/c-sorting-and-keeping-track-of-indexes
73 std::vector<size_t> sort_indexes(const std::vector<T> &v) {
75 // initialize original index locations
76 std::vector<size_t> idx(v.size());
77 iota(idx.begin(), idx.end(), 0);
79 // sort indexes based on comparing values in v
80 std::sort(idx.begin(), idx.end(),
81 [&v](size_t i1, size_t i2) {return v[i1] > v[i2];});
85 //--------------------------------------------------------------------
88 //--------------------------------------------------------------------
89 // Update with the number of dimensions and the pixeltype
90 //--------------------------------------------------------------------
91 template<class args_info_type>
92 template<class InputImageType>
94 MaskOfIntegratedIntensityGenericFilter<args_info_type>::UpdateWithInputImageType()
97 typedef typename InputImageType::PixelType InputPixelType;
98 typedef itk::Image<unsigned char, InputImageType::ImageDimension> MaskImageType;
101 typename InputImageType::Pointer input = this->template GetInput<InputImageType>(0);
103 typename MaskImageType::Pointer mask;
104 mask = MaskImageType::New();
105 mask->SetRegions(input->GetLargestPossibleRegion());
106 mask->SetOrigin(input->GetOrigin());
107 mask->SetSpacing(input->GetSpacing());
111 // Get a vector of all values (will be easier to sort)
112 // And compute total sum of values
113 std::vector<double> values;
114 typedef itk::ImageRegionIterator<InputImageType> IteratorInputType;
115 IteratorInputType iter(input, input->GetLargestPossibleRegion());
118 while (!iter.IsAtEnd()) {
119 values.push_back(iter.Get());
125 auto indices = sort_indexes(values);
127 // Get max index of pixel to reach xx percent
128 double current = 0.0;
129 double max = GetPercentage()/100.0*total;
131 auto n = input->GetLargestPossibleRegion().GetNumberOfPixels();
132 std::vector<int> should_keep(values.size());;
133 std::fill(should_keep.begin(), should_keep.end(), 0);
134 while (current<max and i<n) { // loop by decreasing pixel values
135 current += values[indices[i]];
136 should_keep[indices[i]] = 1.0;
142 typedef itk::ImageRegionIterator<MaskImageType> IteratorMaskType;
143 IteratorMaskType itm(mask, mask->GetLargestPossibleRegion());
147 while (!iter.IsAtEnd()) {
148 if (should_keep[i]) itm.Set(1);
155 if (this->m_IOVerbose)
156 std::cout << "Sum of pixel values : " << total << std::endl
157 << "Percentage : " << GetPercentage() << "%" << std::endl
158 << "Number of pixels : " << nb << "/" << n << std::endl
159 << "Number of pixels : " << nb/n*100.0 << "%" << std::endl;
161 // Write/Save results
162 this->template SetNextOutput<MaskImageType>(mask);
164 //--------------------------------------------------------------------
169 #endif //#define clitkMaskOfIntegratedIntensityGenericFilter_txx