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():
23 // clitk::FilterBase(),
24 clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
25 itk::ImageToImageFilter<ImageType, ImageType>()
27 this->SetNumberOfRequiredInputs(3); // support, object, target
29 SetBackgroundValue(0);
30 SetForegroundValue(1);
33 SetAreaLossTolerance(0.01);
34 m_ListOfAngles.clear();
37 SetSizeWithThreshold(0);
38 SetSizeWithReverseThreshold(0);
40 //--------------------------------------------------------------------
43 //--------------------------------------------------------------------
44 template <class ImageType>
46 clitk::RelativePositionAnalyzerFilter<ImageType>::
47 SetInputSupport(const ImageType * image)
49 // Process object is not const-correct so the const casting is required.
50 this->SetNthInput(0, const_cast<ImageType *>(image));
52 //--------------------------------------------------------------------
55 //--------------------------------------------------------------------
56 template <class ImageType>
58 clitk::RelativePositionAnalyzerFilter<ImageType>::
59 SetInputObject(const ImageType * image)
61 // Process object is not const-correct so the const casting is required.
62 this->SetNthInput(1, const_cast<ImageType *>(image));
64 //--------------------------------------------------------------------
67 //--------------------------------------------------------------------
68 template <class ImageType>
70 clitk::RelativePositionAnalyzerFilter<ImageType>::
71 SetInputTarget(const ImageType * image)
73 // Process object is not const-correct so the const casting is required.
74 this->SetNthInput(2, const_cast<ImageType *>(image));
76 //--------------------------------------------------------------------
79 //--------------------------------------------------------------------
80 template <class ImageType>
82 clitk::RelativePositionAnalyzerFilter<ImageType>::
87 //--------------------------------------------------------------------
90 //--------------------------------------------------------------------
91 template <class ImageType>
93 clitk::RelativePositionAnalyzerFilter<ImageType>::
94 GenerateOutputInformation()
96 ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
97 ImagePointer outputImage = this->GetOutput(0);
98 outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
100 //--------------------------------------------------------------------
103 //--------------------------------------------------------------------
104 template <class ImageType>
106 clitk::RelativePositionAnalyzerFilter<ImageType>::
112 m_Support = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
113 m_Object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
114 m_Target = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(2));
115 static const unsigned int dim = ImageType::ImageDimension;
117 // Remove object from support
118 clitk::AndNot<ImageType>(m_Support, m_Object, GetBackgroundValue());
120 // Resize object like target (to enable substraction later)
121 ImagePointer objectLikeTarget = clitk::ResizeImageLike<ImageType>(m_Object, m_Target, GetBackgroundValue());
123 // Define filter to compute statics on mask image
124 typedef itk::LabelStatisticsImageFilter<ImageType, ImageType> StatFilterType;
125 typename StatFilterType::Pointer statFilter = StatFilterType::New();
127 // Compute the initial support size
128 statFilter->SetInput(m_Support);
129 statFilter->SetLabelInput(m_Support);
130 statFilter->Update();
131 SetSupportSize(statFilter->GetCount(GetForegroundValue()));
132 // DD(GetSupportSize());
134 // Compute the initial target size
135 ImagePointer s = clitk::ResizeImageLike<ImageType>(m_Support, m_Target, GetBackgroundValue());
136 statFilter->SetInput(s);
137 statFilter->SetLabelInput(m_Target);
138 statFilter->Update();
139 SetTargetSize(statFilter->GetCount(GetForegroundValue()));
140 // DD(GetTargetSize());
142 // Build the list of tested orientations
143 m_ListOfAngles.clear();
144 for(uint i=0; i<GetNumberOfAngles(); i++) {
145 double a = i*360.0/GetNumberOfAngles();
146 if (a>180) a = 180-a;
147 m_ListOfAngles.push_back(clitk::deg2rad(a));
148 RelativePositionOrientationType r;
149 r.angle1 = clitk::deg2rad(a);
152 m_ListOfOrientation.push_back(r);
154 m_ListOfOrientation.push_back(r);
157 // Loop on all orientations
158 int bins = GetNumberOfBins();
159 double tolerance = GetAreaLossTolerance();
160 for(int i=0; i<m_ListOfAngles.size(); i++) {
162 typename FloatImageType::Pointer map = ComputeFuzzyMap(objectLikeTarget, m_Target, m_ListOfAngles[i]);
163 writeImage<FloatImageType>(map, "fuzzy_"+toString(i)+".mha");
165 // Compute the optimal thresholds (direct and inverse)
166 double mThreshold=0.0;
167 double mReverseThreshold=1.0;
168 ComputeOptimalThresholds(map, m_Target, bins, tolerance, mThreshold, mReverseThreshold);
170 // Use the threshold to compute new support
171 int s1 = GetSupportSize();
173 // DD(mReverseThreshold);
174 if (mThreshold > 0.0) {
175 ImagePointer support1 =
176 clitk::SliceBySliceRelativePosition<ImageType>(m_Support, m_Object, 2,
178 m_ListOfAngles[i],false,
179 false, -1, true, false);
180 // writeImage<ImageType>(support1, "sup_"+toString(i)+".mha");
181 // Compute the new support size
182 statFilter->SetInput(support1);
183 statFilter->SetLabelInput(support1);
184 statFilter->Update();
185 s1 = statFilter->GetCount(GetForegroundValue());
188 int s2 = GetSupportSize();
189 if (mReverseThreshold < 1.0) {
190 // DD(m_ListOfAngles[1]);
191 ImagePointer support2 =
192 clitk::SliceBySliceRelativePosition<ImageType>(m_Support, m_Object, 2,
194 m_ListOfAngles[i],true,
195 false, -1, true, false);
196 writeImage<ImageType>(support2, "sup_rev_"+toString(i)+".mha");
197 // Compute the new support size
198 statFilter = StatFilterType::New();
199 statFilter->SetInput(support2);
200 statFilter->SetLabelInput(support2);
201 statFilter->Update();
202 s2 = statFilter->GetCount(GetForegroundValue());
205 // Set results values
206 RelativePositionInformationType r;
207 r.threshold = mThreshold;
208 r.sizeAfterThreshold = s1; // DD(s1);
209 r.sizeBeforeThreshold = GetSupportSize();
210 r.sizeReference = GetTargetSize();
211 m_ListOfInformation.push_back(r);
213 r.threshold = mReverseThreshold;
214 r.sizeAfterThreshold = s2; // DD(s2);
215 m_ListOfInformation.push_back(r);
217 } // end loop on orientations
219 //--------------------------------------------------------------------
222 //--------------------------------------------------------------------
223 template <class ImageType>
224 typename clitk::RelativePositionAnalyzerFilter<ImageType>::FloatImageType::Pointer
225 clitk::RelativePositionAnalyzerFilter<ImageType>::
226 ComputeFuzzyMap(ImageType * object, ImageType * target, double angle)
228 typedef clitk::SliceBySliceRelativePositionFilter<ImageType> SliceRelPosFilterType;
229 typedef typename SliceRelPosFilterType::FloatImageType FloatImageType;
230 typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
231 sliceRelPosFilter->VerboseStepFlagOff();
232 sliceRelPosFilter->WriteStepFlagOff();
233 sliceRelPosFilter->SetInput(target);
234 sliceRelPosFilter->SetInputObject(object);
235 sliceRelPosFilter->SetDirection(2);
236 sliceRelPosFilter->SetIntermediateSpacingFlag(false);
237 //sliceRelPosFilter->AddOrientationTypeString(orientation);
238 sliceRelPosFilter->AddAngles(angle, 0.0);
239 sliceRelPosFilter->FuzzyMapOnlyFlagOn(); // do not threshold, only compute the fuzzy map
240 // sliceRelPosFilter->PrintOptions();
241 sliceRelPosFilter->Update();
242 typename FloatImageType::Pointer map = sliceRelPosFilter->GetFuzzyMap();
244 // Remove initial object from the fuzzy map
245 map = clitk::SetBackground<FloatImageType, ImageType>(map, object, GetForegroundValue(), 0.0, true);
247 // Resize the fuzzy map like the target, put 2.0 when outside
248 map = clitk::ResizeImageLike<FloatImageType>(map, target, 2.0); // Put 2.0 when out of initial map
253 //--------------------------------------------------------------------
256 //--------------------------------------------------------------------
257 template <class ImageType>
259 clitk::RelativePositionAnalyzerFilter<ImageType>::
260 ComputeOptimalThresholds(FloatImageType * map, ImageType * target, int bins, double tolerance,
261 double & threshold, double & reverseThreshold)
263 // Get the histogram of fuzzy values inside the target image
264 typedef itk::LabelStatisticsImageFilter<FloatImageType, ImageType> FloatStatFilterType;
265 typename FloatStatFilterType::Pointer f = FloatStatFilterType::New();
267 f->SetLabelInput(target);
268 f->UseHistogramsOn();
269 f->SetHistogramParameters(bins, 0.0, 1.1);
271 int count = f->GetCount(GetForegroundValue());
272 typename FloatStatFilterType::HistogramPointer h = f->GetHistogram(GetForegroundValue());
274 // Debug : dump histogram
276 std::ofstream histogramFile(std::string("fuzzy_histo_"+toString(i)+".txt").c_str());
277 for(int j=0; j<bins; j++) {
278 histogramFile << h->GetMeasurement(j,0)
279 << "\t" << h->GetFrequency(j)
280 << "\t" << (double)h->GetFrequency(j)/(double)count << std::endl;
282 histogramFile.close();
285 // Analyze the histogram (direct)
289 for(int j=0; j<bins; j++) {
290 sum += ((double)h->GetFrequency(j)/(double)count);
291 if ((!found) && (sum > tolerance)) {
292 if (j==0) threshold = h->GetBinMin(0,j);
293 else threshold = h->GetBinMin(0,j-1); // the last before reaching the threshold
298 // Analyze the histogram (reverse)
301 reverseThreshold = 1.0;
302 for(int j=bins-1; j>=0; j--) {
303 sum += ((double)h->GetFrequency(j)/(double)count);
304 if ((!found) && (sum > tolerance)) {
305 if (j==bins-1) reverseThreshold = h->GetBinMax(0,j);
306 else reverseThreshold = h->GetBinMax(0,j+1);
311 //--------------------------------------------------------------------
314 //--------------------------------------------------------------------
315 template <class ImageType>
317 clitk::RelativePositionAnalyzerFilter<ImageType>::
318 Print(std::string s, std::ostream & os)
320 for(int i=0; i<m_ListOfOrientation.size(); i++) {
322 m_ListOfOrientation[i].Print(os);
323 m_ListOfInformation[i].Println(os);
326 //--------------------------------------------------------------------