Orfeo Toolbox  3.16
itkVoronoiSegmentationRGBImageFilter.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3 Program: Insight Segmentation & Registration Toolkit
4 Module: $RCSfile: itkVoronoiSegmentationRGBImageFilter.txx,v $
5 Language: C++
6 Date: $Date: 2010-07-07 12:21:43 $
7 Version: $Revision: 1.37 $
8 
9 Copyright (c) Insight Software Consortium. All rights reserved.
10 See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
11 
12 This software is distributed WITHOUT ANY WARRANTY; without even
13 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14 PURPOSE. See the above copyright notices for more information.
15 
16 =========================================================================*/
17 #ifndef __itkVoronoiSegmentationRGBImageFilter_txx
18 #define __itkVoronoiSegmentationRGBImageFilter_txx
19 #include "itkVoronoiSegmentationRGBImageFilter.h"
20 
21 #include "itkImageRegionIteratorWithIndex.h"
22 #include "itkImageRegionConstIteratorWithIndex.h"
23 #include <math.h>
24 
25 namespace itk
26 {
27 
28 /* Constructor: setting of the default values for the parameters. */
29 template <class TInputImage, class TOutputImage>
30 VoronoiSegmentationRGBImageFilter <TInputImage,TOutputImage>::
31 VoronoiSegmentationRGBImageFilter(){
32  unsigned int i;
33  for(i=0;i<6;i++)
34  {
35  m_Mean[i] = 0;
36  m_STD[i] = 0;
37  m_MeanTolerance[i] = 10;
38  m_STDTolerance[i] = 10;
39  m_MeanPercentError[i] = 0.10;
40  m_STDPercentError[i] = 1.5;
41  }
42 //testing RGB for both mean and STD. (default).
43  m_TestMean[0] = 0;
44  m_TestMean[1] = 1;
45  m_TestMean[2] = 2;
46  m_TestSTD[0] = 0;
47  m_TestSTD[1] = 1;
48  m_TestSTD[2] = 2;
49  m_MaxValueOfRGB = 256;
50 }
51 
52 /* Destructor. */
53 template <class TInputImage, class TOutputImage>
54 VoronoiSegmentationRGBImageFilter <TInputImage,TOutputImage>::
55 ~VoronoiSegmentationRGBImageFilter()
56 {
57 }
58 
59 template <class TInputImage, class TOutputImage>
60 void
61 VoronoiSegmentationRGBImageFilter <TInputImage,TOutputImage>::
62 SetMeanPercentError(double x[6])
63 {
64  for(unsigned int i=0;i<6;i++)
65  {
66  m_MeanPercentError[i] = x[i];
67  m_MeanTolerance[i] = vcl_fabs(x[i]*m_Mean[i]);
68  }
69 }
70 
71 template <class TInputImage, class TOutputImage>
72 void
73 VoronoiSegmentationRGBImageFilter <TInputImage,TOutputImage>::
74 SetSTDPercentError(double x[6])
75 {
76  for(unsigned int i=0;i<6;i++)
77  {
78  m_STDPercentError[i] = x[i];
79  m_STDTolerance[i] = x[i]*m_STD[i];
80  }
81 }
82 
83 
84 /* Initialization for the segmentation. */
85 template <class TInputImage, class TOutputImage>
86 void
87 VoronoiSegmentationRGBImageFilter <TInputImage,TOutputImage>
88 ::SetInput(unsigned int inputNumber, const InputImageType *input)
89 {
90  this->Superclass::SetInput(inputNumber,input);
91 }
92 
93 
94 /* Initialization for the segmentation. */
95 template <class TInputImage, class TOutputImage>
96 void
97 VoronoiSegmentationRGBImageFilter <TInputImage,TOutputImage>
98 ::SetInput(const InputImageType *input)
99 {
100 
101  this->Superclass::SetInput(input);
102 
103  this->SetSize(this->GetInput()->GetLargestPossibleRegion().GetSize());
104  IndexType index;
105  index.Fill(0);
106  RegionType region;
107  region.SetSize(this->GetSize());
108  region.SetIndex(index);
109 
110  m_WorkingImage=RGBHCVImage::New();
111  m_WorkingImage->SetLargestPossibleRegion( region );
112  m_WorkingImage->SetBufferedRegion( region );
113  m_WorkingImage->SetRequestedRegion( region );
114  m_WorkingImage->Allocate();
115 
116  itk::ImageRegionIteratorWithIndex <RGBHCVImage> wit(m_WorkingImage, region);
117  itk::ImageRegionConstIteratorWithIndex <InputImageType> iit(this->GetInput(), region);
118  PixelType ipixel;
119  RGBHCVPixel wpixel;
120 
121  double X;
122  double Y;
123  double Z;
124  double L;
125  double a;
126  double b;
127  double X0 = m_MaxValueOfRGB*0.982;
128  double Y0 = m_MaxValueOfRGB;
129  double Z0 = m_MaxValueOfRGB*1.183;
130 
131  while( !iit.IsAtEnd())
132  {
133  ipixel = iit.Get();
134  wpixel[0] = ipixel[0];
135  wpixel[1] = ipixel[1];
136  wpixel[2] = ipixel[2];
137 
138  X = 0.607*ipixel[0] + 0.174*ipixel[1] + 0.201*ipixel[2];
139  Y = 0.299*ipixel[0] + 0.587*ipixel[1] + 0.114*ipixel[2];
140  Z = 0.066*ipixel[1] + 1.117*ipixel[2];
141  X = vcl_pow((X/X0),0.3333);
142  Y = vcl_pow((Y/Y0),0.3333);
143  Z = vcl_pow((Z/Z0),0.3333);
144  L = 116*Y - 16;
145  a = 500*(X - Y);
146  b = 200*(Y - Z);
147 
148  if (b != 0.0)
149  {
150  wpixel[3] = vcl_atan2(b,a); //H
151  }
152  else
153  {
154  wpixel[3] = 0;
155  }
156  wpixel[4] = vcl_sqrt(a*a+b*b); //C
157  wpixel[5] = L; //V
158  wit.Set(wpixel);
159  ++wit;
160  ++iit;
161  }
162 }
163 
164 template <class TInputImage, class TOutputImage>
165 bool
166 VoronoiSegmentationRGBImageFilter <TInputImage,TOutputImage>::
167 TestHomogeneity(IndexList &Plist)
168 {
169  int num=Plist.size();
170  int i,j;
171  RGBHCVPixel getp;
172  double addp[6]={0,0,0,0,0,0};
173  double addpp[6]={0,0,0,0,0,0};
174  for(i=0;i<num;i++)
175  {
176  getp = m_WorkingImage->GetPixel(Plist[i]);
177  for(j=0;j<6;j++)
178  {
179  addp[j]=addp[j]+getp[j];
180  addpp[j]=addpp[j]+getp[j]*getp[j];
181  }
182  }
183 
184  double savemean[6],saveSTD[6];
185  if(num > 1)
186  {
187  for(i=0;i<6;i++)
188  {
189  savemean[i] = addp[i]/num;
190  saveSTD[i] = vcl_sqrt((addpp[i] - (addp[i]*addp[i])/(num) )/(num-1));
191  }
192  }
193  else
194  {
195  for(i=0;i<6;i++)
196  {
197  savemean[i] = 0;
198  saveSTD[i] = -1;
199  }
200  }
201 
202 
203  bool ok = 1;
204  j = 0;
205  double savem,savev;
206  while (ok && (j < 3))
207  {
208  savem = savemean[m_TestMean[j]] - m_Mean[m_TestMean[j]];
209  savev = saveSTD[m_TestSTD[j]] - m_STD[m_TestSTD[j]];
210  if( (savem < -m_MeanTolerance[m_TestMean[j]]) ||
211  (savem > m_MeanTolerance[m_TestMean[j]]) )
212  {
213  ok = 0;
214  }
215  if( (savev < -m_STDTolerance[m_TestSTD[j]]) ||
216  (savev > m_STDTolerance[m_TestSTD[j]]) )
217  {
218  ok = 0;
219  }
220  j++;
221  }
222 
223  if(ok)
224  return 1;
225  else
226  return 0;
227 }
228 
229 template <class TInputImage, class TOutputImage>
230 void
231 VoronoiSegmentationRGBImageFilter <TInputImage,TOutputImage>::
232 TakeAPrior(const BinaryObjectImage* aprior)
233 {
234 
235  RegionType region = this->GetInput()->GetRequestedRegion();
236  itk::ImageRegionConstIteratorWithIndex <BinaryObjectImage> ait(aprior, region);
237  itk::ImageRegionIteratorWithIndex <RGBHCVImage> iit(m_WorkingImage, region);
238 
239  unsigned int minx=0,miny=0,maxx=0,maxy=0;
240  bool status=0;
241  for(unsigned int i=0;i<this->GetSize()[1];i++)
242  {
243  for(unsigned int j=0;j<this->GetSize()[0];j++)
244  {
245  if( (status==0)&&(ait.Get()) )
246  {
247  miny=i;
248  minx=j;
249  maxy=i;
250  maxx=j;
251  status=1;
252  }
253  else if( (status==1)&&(ait.Get()) )
254  {
255  maxy=i;
256  if(minx>j) minx=j;
257  if(maxx<j) maxx=j;
258  }
259  ++ait;
260  }
261  }
262 
263  int objnum = 0;
264  int bkgnum = 0;
265 
266  float objaddp[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
267  float objaddpp[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
268  float bkgaddp[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
269  float bkgaddpp[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
270  RGBHCVPixel currp;
271 
272  ait.GoToBegin();
273  iit.GoToBegin();
274  unsigned int k;
275  for(unsigned int i=0;i<miny;i++)
276  {
277  for(unsigned int j=0;j<this->GetSize()[0];j++)
278  {
279  ++ait;
280  ++iit;
281  }
282  }
283  for(unsigned int i=miny;i<=maxy;i++)
284  {
285  for(unsigned int j=0;j<minx;j++)
286  {
287  ++ait;
288  ++iit;
289  }
290  for(unsigned int j=minx;j<=maxx;j++)
291  {
292  currp = iit.Get();
293  if(ait.Get())
294  {
295  objnum++;
296  for(k=0;k<6;k++)
297  {
298  objaddp[k] += currp[k];
299  objaddpp[k] += currp[k]*currp[k];
300  }
301  }
302  else
303  {
304  bkgnum++;
305  for(k=0;k<6;k++)
306  {
307  bkgaddp[k] += currp[k];
308  bkgaddpp[k] += currp[k]*currp[k];
309  }
310  }
311  ++ait;++iit;
312  }
313  for(unsigned int j=maxx+1;j<this->GetSize()[0];j++)
314  {
315  ++ait;
316  ++iit;
317  }
318  }
319 
320  double b_Mean[6];
321  double b_STD[6];
322  float diffMean[6];
323  float diffSTD[6];
324  for(unsigned int i=0;i<6;i++)
325  {
326  m_Mean[i] = objaddp[i]/objnum;
327  m_STD[i] = vcl_sqrt((objaddpp[i] - (objaddp[i]*objaddp[i])/objnum)/(objnum-1));
328  m_STDTolerance[i] = m_STD[i]*m_STDPercentError[i];
329  b_Mean[i] = bkgaddp[i]/bkgnum;
330  b_STD[i] = vcl_sqrt((bkgaddpp[i] - (bkgaddp[i]*bkgaddp[i])/bkgnum)/(bkgnum-1));
331  diffMean[i] = (b_Mean[i]-m_Mean[i])/m_Mean[i];
332  if(diffMean[i] < 0) diffMean[i] = -diffMean[i];
333  diffSTD[i] = (b_STD[i]-m_STD[i])/m_STD[i];
334  if(diffSTD[i] < 0) diffSTD[i] = -diffSTD[i];
335  if(this->GetUseBackgroundInAPrior())
336  {
337  m_MeanTolerance[i] = diffMean[i]*m_Mean[i]*this->GetMeanDeviation();
338  }
339  else
340  {
341  m_MeanTolerance[i] = vcl_fabs(m_Mean[i]*m_MeanPercentError[i]);
342  }
343  }
344 
345  if(objnum<10)
346  {
347 /* a-prior doen's make too much sense */
348  for(unsigned int i=0;i<6;i++)
349  {
350  m_MeanTolerance[i] = 0;
351  m_STDTolerance[i] = 0;
352  }
353  }
354 
355 /* Sorting. */
356  unsigned char tmp[6]={0,1,2,3,4,5};
357  for(unsigned j=0;j<3;j++)
358  {
359  k=0;
360  for(unsigned int i=1;i<6-j;i++)
361  {
362  if(diffMean[tmp[i]]>diffMean[tmp[k]])
363  {
364  k=i;
365  }
366  }
367  m_TestMean[j]=tmp[k];
368  tmp[k]=tmp[5-j];
369  }
370  unsigned char tmp1[6]={0,1,2,3,4,5};
371  for(unsigned int j=0;j<3;j++)
372  {
373  k=0;
374  for(unsigned int i=1;i<6-j;i++)
375  {
376  if(diffSTD[tmp1[i]]>diffSTD[tmp1[k]])
377  {
378  k=i;
379  }
380  }
381  m_TestSTD[j]=tmp1[k];
382  tmp1[k]=tmp1[5-j];
383  }
384 }
385 template <class TInputImage, class TOutputImage>
386 void
387 VoronoiSegmentationRGBImageFilter <TInputImage,TOutputImage>
388 ::PrintSelf(std::ostream& os, Indent indent) const
389 {
390  Superclass::PrintSelf(os, indent);
391  os << indent << "MaxValueOfRGB: " << m_MaxValueOfRGB << std::endl;
392  os << indent << "Mean: " << m_Mean << std::endl;
393 
394 }
395 } //end namespace
396 
397 #endif
Generated at Sun Feb 3 2013 00:13:55 for Orfeo Toolbox with doxygen 1.8.1.1