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 //--------------------------------------------------------------------
20 template <class ImageType>
21 clitk::RelativePositionAnalyzerFilter<ImageType>::
22 RelativePositionAnalyzerFilter():
24 clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
25 itk::ImageToImageFilter<ImageType, ImageType>()
27 this->SetNumberOfRequiredInputs(3); // support, object, target
29 SetBackgroundValue(0);
30 SetForegroundValue(1);
32 //--------------------------------------------------------------------
35 //--------------------------------------------------------------------
36 template <class ImageType>
38 clitk::RelativePositionAnalyzerFilter<ImageType>::
39 SetInputSupport(const ImageType * image)
41 // Process object is not const-correct so the const casting is required.
42 this->SetNthInput(0, const_cast<ImageType *>(image));
44 //--------------------------------------------------------------------
47 //--------------------------------------------------------------------
48 template <class ImageType>
50 clitk::RelativePositionAnalyzerFilter<ImageType>::
51 SetInputObject(const ImageType * image)
53 // Process object is not const-correct so the const casting is required.
54 this->SetNthInput(1, const_cast<ImageType *>(image));
56 //--------------------------------------------------------------------
59 //--------------------------------------------------------------------
60 template <class ImageType>
62 clitk::RelativePositionAnalyzerFilter<ImageType>::
63 SetInputTarget(const ImageType * image)
65 // Process object is not const-correct so the const casting is required.
66 this->SetNthInput(2, const_cast<ImageType *>(image));
68 //--------------------------------------------------------------------
71 //--------------------------------------------------------------------
72 template <class ImageType>
74 clitk::RelativePositionAnalyzerFilter<ImageType>::
79 //--------------------------------------------------------------------
82 //--------------------------------------------------------------------
83 template <class ImageType>
85 clitk::RelativePositionAnalyzerFilter<ImageType>::
86 GenerateOutputInformation()
88 ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
89 ImagePointer outputImage = this->GetOutput(0);
90 outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
92 //--------------------------------------------------------------------
95 //--------------------------------------------------------------------
96 template <class ImageType>
98 clitk::RelativePositionAnalyzerFilter<ImageType>::
99 GenerateInputRequestedRegion()
102 itk::ImageToImageFilter<ImageType, ImageType>::GenerateInputRequestedRegion();
103 // Get input pointers and set requested region to common region
104 ImagePointer input1 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
105 ImagePointer input2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
106 ImagePointer input3 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(2));
107 input1->SetRequestedRegion(input1->GetLargestPossibleRegion());
108 input2->SetRequestedRegion(input2->GetLargestPossibleRegion());
109 input3->SetRequestedRegion(input3->GetLargestPossibleRegion());
111 //--------------------------------------------------------------------
113 //--------------------------------------------------------------------
114 template <class ImageType>
116 clitk::RelativePositionAnalyzerFilter<ImageType>::
122 m_Support = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
123 m_Object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
124 m_Target = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(2));
125 static const unsigned int dim = ImageType::ImageDimension;
127 // Remove object from support
128 clitk::AndNot<ImageType>(m_Support, m_Object, GetBackgroundValue());
130 // Resize object like target (to enable substraction later)
131 ImagePointer objectLikeTarget = clitk::ResizeImageLike<ImageType>(m_Object, m_Target, GetBackgroundValue());
133 // Define filter to compute statics on mask image
134 typedef itk::LabelStatisticsImageFilter<ImageType, ImageType> StatFilterType;
135 typename StatFilterType::Pointer statFilter = StatFilterType::New();
137 // Compute the initial support size
138 statFilter->SetInput(m_Support);
139 statFilter->SetLabelInput(m_Support);
140 statFilter->Update();
141 int m_SupportSize = statFilter->GetCount(GetForegroundValue());
142 // DD(m_SupportSize);
144 // Compute the initial target size
145 ImagePointer s = clitk::ResizeImageLike<ImageType>(m_Support, m_Target, GetBackgroundValue());
146 statFilter->SetInput(s);
147 statFilter->SetLabelInput(m_Target);
148 statFilter->Update();
149 int m_TargetSize = statFilter->GetCount(GetForegroundValue());
152 // Build the list of tested orientations
153 std::vector<double> m_ListOfAngles;
154 int m_NumberOfAngles = this->GetAFDB()->GetDouble("NumberOfAngles");
155 for(uint i=0; i<m_NumberOfAngles; i++) {
156 double a = i*360.0/m_NumberOfAngles;
157 if (a>180) a = 180-a;
158 m_ListOfAngles.push_back(clitk::deg2rad(a));
161 // Loop on all orientations
162 for(int i=0; i<m_ListOfAngles.size(); i++) {
165 typename FloatImageType::Pointer map = ComputeFuzzyMap(objectLikeTarget, m_Target, m_ListOfAngles[i]);
166 writeImage<FloatImageType>(map, "fuzzy_"+toString(i)+".mha");
168 // Compute the optimal thresholds (direct and inverse)
170 double mReverseThreshold;
171 int bins = this->GetAFDB()->GetDouble("bins");
172 double tolerance = this->GetAFDB()->GetDouble("TargetAreaLossTolerance");
173 ComputeOptimalThresholds(map, m_Target, bins, tolerance, mThreshold, mReverseThreshold);
175 // Use the threshold to compute new support
177 if (mThreshold > 0.0) {
178 ImagePointer support1 =
179 clitk::SliceBySliceRelativePosition<ImageType>(m_Support, m_Object, 2,
181 m_ListOfAngles[i],false,
182 false, -1, true, false);
183 writeImage<ImageType>(support1, "sup_"+toString(i)+".mha");
184 // Compute the new support size
185 statFilter->SetInput(support1);
186 statFilter->SetLabelInput(support1);
187 statFilter->Update();
188 s1 = statFilter->GetCount(GetForegroundValue());
190 else s1 = m_SupportSize;
193 if (mReverseThreshold < 1.0) {
194 ImagePointer support2 =
195 clitk::SliceBySliceRelativePosition<ImageType>(m_Support, m_Object, 2,
197 m_ListOfAngles[i],true,
198 false, -1, true, false);
199 writeImage<ImageType>(support2, "sup_rev_"+toString(i)+".mha");
200 // Compute the new support size
201 statFilter = StatFilterType::New();
202 statFilter->SetInput(support2);
203 statFilter->SetLabelInput(support2);
204 statFilter->Update();
205 s2 = statFilter->GetCount(GetForegroundValue());
207 else s2 =m_SupportSize;
210 std::cout << i << " " << clitk::rad2deg(m_ListOfAngles[i]) << "\t"
211 << m_SupportSize << " " << m_TargetSize << "\t"
212 << s1/(double)m_SupportSize << " " << s2/(double)m_SupportSize << "\t"
213 << mThreshold << " " << mReverseThreshold << std::endl;
215 } // end loop on orientations
218 // Final Step -> set output TODO
219 // this->SetNthOutput(0, working_image);
220 // this->GraftOutput(working_image);
222 //--------------------------------------------------------------------
225 //--------------------------------------------------------------------
226 template <class ImageType>
227 typename clitk::RelativePositionAnalyzerFilter<ImageType>::FloatImageType::Pointer
228 clitk::RelativePositionAnalyzerFilter<ImageType>::
229 ComputeFuzzyMap(ImageType * object, ImageType * target, double angle)
231 typedef clitk::SliceBySliceRelativePositionFilter<ImageType> SliceRelPosFilterType;
232 typedef typename SliceRelPosFilterType::FloatImageType FloatImageType;
233 typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
234 sliceRelPosFilter->VerboseStepFlagOff();
235 sliceRelPosFilter->WriteStepFlagOff();
236 sliceRelPosFilter->SetInput(target);
237 sliceRelPosFilter->SetInputObject(object);
238 sliceRelPosFilter->SetDirection(2);
239 sliceRelPosFilter->SetIntermediateSpacingFlag(false);
240 //sliceRelPosFilter->AddOrientationTypeString(orientation);
241 sliceRelPosFilter->AddAngles(angle, 0.0);
242 sliceRelPosFilter->FuzzyMapOnlyFlagOn(); // do not threshold, only compute the fuzzy map
243 // sliceRelPosFilter->PrintOptions();
244 sliceRelPosFilter->Update();
245 typename FloatImageType::Pointer map = sliceRelPosFilter->GetFuzzyMap();
247 // Remove initial object from the fuzzy map
248 map = clitk::SetBackground<FloatImageType, ImageType>(map, object, GetForegroundValue(), 0.0, true);
250 // Resize the fuzzy map like the target, put 2.0 when outside
251 map = clitk::ResizeImageLike<FloatImageType>(map, target, 2.0); // Put 2.0 when out of initial map
256 //--------------------------------------------------------------------
259 //--------------------------------------------------------------------
260 template <class ImageType>
262 clitk::RelativePositionAnalyzerFilter<ImageType>::
263 ComputeOptimalThresholds(FloatImageType * map, ImageType * target, int bins, double tolerance,
264 double & threshold, double & reverseThreshold)
266 // Get the histogram of fuzzy values inside the target image
267 typedef itk::LabelStatisticsImageFilter<FloatImageType, ImageType> FloatStatFilterType;
268 typename FloatStatFilterType::Pointer f = FloatStatFilterType::New();
270 f->SetLabelInput(target);
271 f->UseHistogramsOn();
272 f->SetHistogramParameters(bins, 0.0, 1.1);
274 int count = f->GetCount(GetForegroundValue());
275 typename FloatStatFilterType::HistogramPointer h = f->GetHistogram(GetForegroundValue());
277 // Debug : dump histogram
279 std::ofstream histogramFile(std::string("fuzzy_histo_"+toString(i)+".txt").c_str());
280 for(int j=0; j<bins; j++) {
281 histogramFile << h->GetMeasurement(j,0)
282 << "\t" << h->GetFrequency(j)
283 << "\t" << (double)h->GetFrequency(j)/(double)count << std::endl;
285 histogramFile.close();
288 // Analyze the histogram (direct)
292 for(int j=0; j<bins; j++) {
293 sum += ((double)h->GetFrequency(j)/(double)count);
294 if ((!found) && (sum > tolerance)) {
295 threshold = h->GetBinMin(0,j);
300 // Analyze the histogram (reverse)
303 reverseThreshold = 1.0;
304 for(int j=bins-1; j>=0; j--) {
305 sum += ((double)h->GetFrequency(j)/(double)count);
306 if ((!found) && (sum > tolerance)) {
307 reverseThreshold = h->GetBinMax(0,j);
312 //--------------------------------------------------------------------