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 std::vector<std::pair<T, size_t> > compVector(v.size());
78 for (size_t i = 0; i < v.size(); ++i) {
79 compVector[i].first = v[i];
80 compVector[i].second = i;
83 // sort indexes based on comparing values in v
84 std::sort(compVector.begin(), compVector.end(), comparator<T>);
85 for (size_t i = 0; i < v.size(); ++i) {
86 idx[i] = compVector[i].second;
91 //--------------------------------------------------------------------
94 //--------------------------------------------------------------------
95 // Update with the number of dimensions and the pixeltype
96 //--------------------------------------------------------------------
97 template<class args_info_type>
98 template<class InputImageType>
100 MaskOfIntegratedIntensityGenericFilter<args_info_type>::UpdateWithInputImageType()
103 typedef typename InputImageType::PixelType InputPixelType;
104 typedef itk::Image<unsigned char, InputImageType::ImageDimension> MaskImageType;
107 typename InputImageType::Pointer input = this->template GetInput<InputImageType>(0);
109 typename MaskImageType::Pointer mask;
110 mask = MaskImageType::New();
111 mask->SetRegions(input->GetLargestPossibleRegion());
112 mask->SetOrigin(input->GetOrigin());
113 mask->SetSpacing(input->GetSpacing());
117 // Get a vector of all values (will be easier to sort)
118 // And compute total sum of values
119 std::vector<double> values;
120 typedef itk::ImageRegionIterator<InputImageType> IteratorInputType;
121 IteratorInputType iter(input, input->GetLargestPossibleRegion());
124 while (!iter.IsAtEnd()) {
125 values.push_back(iter.Get());
131 std::vector<size_t> indices = sort_indexes(values);
133 // Get max index of pixel to reach xx percent
134 double current = 0.0;
135 double max = GetPercentage()/100.0*total;
137 int n = input->GetLargestPossibleRegion().GetNumberOfPixels();
138 std::vector<int> should_keep(values.size());;
139 std::fill(should_keep.begin(), should_keep.end(), 0);
140 while (current<max && i<n) { // loop by decreasing pixel values
141 current += values[indices[i]];
142 should_keep[indices[i]] = 1.0;
148 typedef itk::ImageRegionIterator<MaskImageType> IteratorMaskType;
149 IteratorMaskType itm(mask, mask->GetLargestPossibleRegion());
153 while (!iter.IsAtEnd()) {
154 if (should_keep[i]) itm.Set(1);
161 if (this->m_IOVerbose)
162 std::cout << "Sum of pixel values : " << total << std::endl
163 << "Percentage : " << GetPercentage() << "%" << std::endl
164 << "Number of pixels : " << nb << "/" << n << std::endl
165 << "Number of pixels : " << nb/n*100.0 << "%" << std::endl;
167 // Write/Save results
168 this->template SetNextOutput<MaskImageType>(mask);
170 //--------------------------------------------------------------------
172 //--------------------------------------------------------------------
173 template <typename T>
174 bool comparator ( const std::pair<T, size_t>& l, const std::pair<T, size_t>& r)
175 { return l.first > r.first; }
176 //--------------------------------------------------------------------
180 #endif //#define clitkMaskOfIntegratedIntensityGenericFilter_txx