Orfeo Toolbox  3.16
itkVoronoiSegmentationImageFilterBase.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: Insight Segmentation & Registration Toolkit
4  Module: $RCSfile: itkVoronoiSegmentationImageFilterBase.txx,v $
5  Language: C++
6  Date: $Date: 2009-04-06 13:46:34 $
7  Version: $Revision: 1.35 $
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 __itkVoronoiSegmentationImageFilterBase_txx
18 #define __itkVoronoiSegmentationImageFilterBase_txx
19 
20 #include "itkImageRegionIteratorWithIndex.h"
21 #include "itkVoronoiDiagram2DGenerator.h"
22 #include <math.h>
23 
24 namespace itk
25 {
26 
27 /* Constructor: setting the default parameter values. */
28 template <class TInputImage, class TOutputImage, class TBinaryPriorImage>
29 VoronoiSegmentationImageFilterBase <TInputImage,TOutputImage,TBinaryPriorImage>
30 ::VoronoiSegmentationImageFilterBase()
31 {
32  m_Size.Fill(0);
33  m_MinRegion = 20;
34  m_Steps = 0;
35  m_LastStepSeeds = 0;
36  m_NumberOfSeeds = 200;
37  m_NumberOfSeedsToAdded = 0;
38  m_MeanDeviation = 0.8;
39  m_UseBackgroundInAPrior = false;
40  m_OutputBoundary = false;
41  m_InteractiveSegmentation = false;
42  m_WorkingVD=VoronoiDiagram::New();
43  m_VDGenerator=VoronoiDiagramGenerator::New();
44 }
45 
46 /* Destructor. */
47 template <class TInputImage, class TOutputImage, class TBinaryPriorImage>
48 VoronoiSegmentationImageFilterBase <TInputImage,TOutputImage,TBinaryPriorImage>
49 ::~VoronoiSegmentationImageFilterBase()
50 {
51 }
52 
53 template <class TInputImage, class TOutputImage, class TBinaryPriorImage>
54 void
55 VoronoiSegmentationImageFilterBase <TInputImage,TOutputImage,TBinaryPriorImage>
56 ::PrintSelf(std::ostream& os, Indent indent) const
57 {
58  Superclass::PrintSelf(os, indent);
59  os << indent << "Number Of Seeds: "
60  << m_NumberOfSeeds << std::endl;
61 
62  os << indent << "Minimum Region for Split: "
63  << m_MinRegion << std::endl;
64 
65  os << indent << "Number Of Steps to Run: (0 means runs until no region to split) "
66  << m_Steps << std::endl;
67 
68  os << indent << "UseBackgroundInAPrior = " << m_UseBackgroundInAPrior << std::endl;
69  os << indent << "OutputBoundary = " << m_OutputBoundary << std::endl;
70  os << indent << "MeanDeviation = " << m_MeanDeviation << std::endl;
71  os << indent << "LastStepSeeds = " << m_LastStepSeeds << std::endl;
72  os << indent << "InteractiveSegmentation = " << m_InteractiveSegmentation << std::endl;
73  os << indent << "NumberOfSeedsToAdded = " << m_NumberOfSeedsToAdded
74  << std::endl;
75  os << indent << "Size = " << m_Size << std::endl;
76 }
77 
78 template <class TInputImage, class TOutputImage, class TBinaryPriorImage>
79 void
80 VoronoiSegmentationImageFilterBase <TInputImage,TOutputImage,TBinaryPriorImage>
81 ::GetPixelIndexFromPolygon(PointTypeDeque vertlist, IndexList *PixelPool)
82 {
83  IndexType idx;
84  PointType currP;
85  PointType leftP;
86  PointType rightP;
87 
88  currP = vertlist.front();
89  vertlist.pop_front();
90  leftP = vertlist.front();
91  while(currP[1] > leftP[1])
92  {
93  vertlist.push_back(currP);
94  currP=vertlist.front();
95  vertlist.pop_front();
96  leftP=vertlist.front();
97  }
98  rightP=vertlist.back();
99  while(currP[1] > rightP[1])
100  {
101  vertlist.push_front(currP);
102  currP=vertlist.back();
103  vertlist.pop_back();
104  rightP=vertlist.back();
105  }
106  leftP=vertlist.front();
107  PointTypeDeque tmpQ;
108  tmpQ.clear();
109  if(leftP[0]>rightP[0])
110  {
111  while(!(vertlist.empty()))
112  {
113  tmpQ.push_back(vertlist.front());
114  vertlist.pop_front();
115  }
116  while(!(tmpQ.empty()))
117  {
118  vertlist.push_front(tmpQ.front());
119  tmpQ.pop_front();
120  }
121  }
122  tmpQ.clear();
123  leftP=vertlist.front();
124  rightP=vertlist.back();
125 
126  double beginy=currP[1];
127  int intbeginy=(int)vcl_ceil(beginy);
128  idx[1]=intbeginy;
129  double leftendy=leftP[1];
130  double rightendy=rightP[1];
131  double beginx=currP[0];
132  double endx=currP[0];
133  double leftDx,rightDx;
134  double offset;
135  double leftheadx=beginx;
136  double rightheadx=endx;
137  double leftheady=beginy;
138  double rightheady=beginy;
139 
140  double endy;
141  bool RorL;
142  int i;
143  if(leftendy>rightendy)
144  {
145  RorL=1;
146  endy=rightendy;
147  }
148  else
149  {
150  RorL=0;
151  endy=leftendy;
152  }
153  if (leftP[1]==beginy)
154  {
155  leftDx = 1.0;
156  }
157  else
158  {
159  leftDx=(leftP[0]-beginx)/(leftP[1]-beginy);
160  }
161 
162  if (rightP[1]==beginy)
163  {
164  rightDx = 1.0;
165  }
166  else
167  {
168  rightDx=(rightP[0]-endx)/(rightP[1]-beginy);
169  }
170 
171  int intendy=(int)vcl_floor(endy);
172  if(intbeginy>intendy)
173  { //no scanline
174  if(RorL)
175  {
176  beginy=rightP[1];
177  }
178  else
179  {
180  beginy=leftP[1];
181  }
182  }
183  else if((intbeginy==intendy) && (intbeginy==0))
184  { //only one scanline at 0;
185  if(RorL)
186  {
187  endx=rightP[0];
188  }
189  else
190  {
191  beginx=leftP[0];
192  }
193  for(i=static_cast<int>(vcl_ceil(beginx));i<=static_cast<int>(vcl_floor(endx));i++)
194  {
195  idx[0]=i;
196  (*PixelPool).push_back(idx);
197  }
198  idx[1]=idx[1]+1;
199  }
200  else
201  { //normal case some scanlines
202  offset = (double)intbeginy-beginy;
203  endx += offset * rightDx;
204  beginx += offset * leftDx;
205  while(idx[1] <= intendy)
206  {
207  for(i=static_cast<int>(vcl_ceil(beginx)); i <= static_cast<int>(vcl_floor(endx)); i++)
208  {
209  idx[0] = i;
210  (*PixelPool).push_back(idx);
211  }
212  endx += rightDx;
213  beginx += leftDx;
214  idx[1] = idx[1]+1;
215  }
216  beginy=endy;
217  }
218 
219  int vsize = static_cast<int>( vertlist.size() );
220  while(vsize>2)
221  {
222  vsize--;
223  if(RorL)
224  {
225  vertlist.pop_back();
226  currP=rightP;
227  rightheadx=currP[0];
228  rightheady=currP[1];
229  endx=currP[0];
230  beginx=leftheadx+leftDx*(beginy-leftheady);
231  rightP=vertlist.back();
232  if ((rightP[1]==currP[1]))
233  {
234  rightDx = 1.0;
235  }
236  else
237  {
238  rightDx=(rightP[0]-currP[0])/(rightP[1]-currP[1]);
239  }
240  }
241  else
242  {
243  vertlist.pop_front();
244  currP=leftP;
245  leftheadx=currP[0];
246  leftheady=currP[1];
247  beginx=currP[0];
248  endx=rightheadx+rightDx*(beginy-rightheady);
249  leftP=vertlist.front();
250  if ((leftP[1]==currP[1]))
251  {
252  leftDx = 1.0;
253  }
254  else
255  {
256  leftDx=(leftP[0]-currP[0])/(leftP[1]-currP[1]);
257  }
258  }
259 
260  leftendy=leftP[1];
261  rightendy=rightP[1];
262  if(leftendy>rightendy)
263  {
264  RorL=1;
265  endy=rightendy;
266  }
267  else
268  {
269  RorL=0;
270  endy=leftendy;
271  }
272 
273  intendy=(int)vcl_floor(endy);
274  intbeginy=(int)vcl_ceil(beginy);
275 
276  if(intbeginy>intendy)
277  { //no scanline
278  if(RorL)
279  {
280  beginy=rightP[1];
281  }
282  else
283  {
284  beginy=leftP[1];
285  }
286  }
287  else
288  { //normal case some scanlines
289  offset = (double)intbeginy-beginy;
290  endx += offset * rightDx;
291  beginx += offset * leftDx;
292  while(idx[1] <= intendy)
293  {
294  for(i = static_cast<int>(vcl_ceil(beginx)); i <= static_cast<int>(vcl_floor(endx)); i++)
295  {
296  idx[0] = i;
297  (*PixelPool).push_back(idx);
298  }
299  endx += rightDx;
300  beginx += leftDx;
301  idx[1] = idx[1] + 1;
302  }
303  beginy=idx[1];
304  }
305  }
306 
307 
308  if(RorL)
309  {
310  beginy=rightP[1];
311  endy=leftP[1];
312  }
313  else
314  {
315  beginy=leftP[1];
316  endy=rightP[1];
317  }
318  intbeginy=(int)vcl_ceil(beginy);
319  intendy=(int)vcl_floor(endy);
320  if(intbeginy<=intendy)
321  {
322  if(RorL)
323  {
324  if ((rightP[1]==leftP[1]))
325  {
326  rightDx = 1.0;
327  }
328  else
329  {
330  rightDx=(rightP[0]-leftP[0])/(rightP[1]-leftP[1]);
331  }
332  endx=rightP[0];
333  beginx=leftP[0]+leftDx*(rightP[1]-leftP[1]);
334  }
335  else
336  {
337  if ((rightP[1]==leftP[1]))
338  {
339  leftDx = 1.0;
340  }
341  else
342  {
343  leftDx=(rightP[0]-leftP[0])/(rightP[1]-leftP[1]);
344  }
345  beginx=leftP[0];
346  endx=rightP[0]+rightDx*(leftP[1]-rightP[1]);
347  }
348  offset = (double)intbeginy-beginy;
349  beginx += offset * leftDx;
350  endx += offset * rightDx;
351  while(idx[1] <= intendy)
352  {
353  for(i = static_cast<int>(vcl_ceil(beginx)); i <= static_cast<int>(vcl_floor(endx)); i++)
354  {
355  idx[0] = i;
356  (*PixelPool).push_back(idx);
357  }
358  endx += rightDx;
359  beginx += leftDx;
360  idx[1] = idx[1] + 1;
361  }
362  }
363 }
364 
365 template <class TInputImage, class TOutputImage, class TBinaryPriorImage>
366 void
367 VoronoiSegmentationImageFilterBase <TInputImage,TOutputImage,TBinaryPriorImage>
368 ::ClassifyDiagram(void)
369 {
370  PointIdIterator currPit;
371  PointIdIterator currPitEnd;
372  PointType currP;
373  PointTypeDeque VertList;
374  IndexList PixelPool;
375  for(int i=0;i<m_NumberOfSeeds;i++)
376  {
377  CellAutoPointer currCell;
378  m_WorkingVD->GetCellId(i,currCell);
379  currPitEnd = currCell->PointIdsEnd();
380  VertList.clear();
381  for(currPit = currCell->PointIdsBegin(); currPit != currPitEnd; ++currPit)
382  {
383  m_WorkingVD->GetPoint((*currPit),&(currP));
384  VertList.push_back(currP);
385  }
386 
387  PixelPool.clear();
388  this->GetPixelIndexFromPolygon(VertList,&PixelPool);
389  m_NumberOfPixels[i] = static_cast<int>( PixelPool.size() );
390  m_Label[i] = this->TestHomogeneity(PixelPool);
391  }
392 
393  m_NumberOfBoundary = 0;
394  for(int i=0;i<m_NumberOfSeeds;i++)
395  {
396  // if not homogeneous
397  if(m_Label[i] == 0)
398  {
399  NeighborIdIterator itend = m_WorkingVD->NeighborIdsEnd(i);
400  NeighborIdIterator it=m_WorkingVD->NeighborIdsBegin(i);
401  bool bnd = 0;
402  // and any adjacent region is homogeneous
403  while((it != itend) && (!bnd))
404  {
405  bnd = (m_Label[*it] == 1);
406  ++it;
407  }
408  // then this must be a boundary region
409  if(bnd)
410  {
411  m_Label[i] = 2;
412  m_NumberOfBoundary++;
413  }
414  }
415  }
416 }
417 
418 template <class TInputImage, class TOutputImage, class TBinaryPriorImage>
419 void
420 VoronoiSegmentationImageFilterBase <TInputImage,TOutputImage,TBinaryPriorImage>
421 ::GenerateAddingSeeds(void)
422 {
423  EdgeIterator eit;
424  EdgeIterator eitend = m_WorkingVD->EdgeEnd();
425  PointType adds;
426  Point<int,2> seeds;
427 
428  for(eit = m_WorkingVD->EdgeBegin();eit != eitend; ++eit)
429  {
430  seeds = m_WorkingVD->GetSeedsIDAroundEdge(&*eit);
431  if( ((m_Label[seeds[0]]==2)||(m_Label[seeds[1]]==2))
432  && (m_NumberOfPixels[seeds[0]]>m_MinRegion)
433  && (m_NumberOfPixels[seeds[1]]>m_MinRegion) )
434  {
435  adds[0] = (eit->m_Left[0] + eit->m_Right[0])*0.5;
436  adds[1] = (eit->m_Left[1] + eit->m_Right[1])*0.5;
437  m_SeedsToAdded.push_back(adds);
438  }
439  }
440 }
441 
442 
443 template <class TInputImage, class TOutputImage, class TBinaryPriorImage>
444 void
445 VoronoiSegmentationImageFilterBase <TInputImage,TOutputImage,TBinaryPriorImage>
446 ::RunSegmentOneStep(void)
447 {
448  m_NumberOfPixels.resize(m_NumberOfSeeds);
449  m_Label.resize(m_NumberOfSeeds);
450  m_SeedsToAdded.clear();
451  m_VDGenerator->Update();
452  m_WorkingVD=m_VDGenerator->GetOutput();
453  this->ClassifyDiagram();
454  this->GenerateAddingSeeds();
455  m_NumberOfSeedsToAdded = static_cast<int>( m_SeedsToAdded.size() );
456 
457  if (m_InteractiveSegmentation)
458  {
459  if(m_OutputBoundary)
460  {
461  this->MakeSegmentBoundary();
462  }
463  else
464  {
465  this->MakeSegmentObject();
466  }
467  }
468 }
469 
470 template <class TInputImage, class TOutputImage, class TBinaryPriorImage>
471 void
472 VoronoiSegmentationImageFilterBase <TInputImage,TOutputImage,TBinaryPriorImage>
473 ::RunSegment(void)
474 {
475  bool ok = 1;
476  if(m_Steps == 0)
477  {
478  unsigned long count=1;
479  this->RunSegmentOneStep();
480  // guess at a progress
481  this->UpdateProgress(static_cast<float>(count)
482  / static_cast<float>(NumericTraits<unsigned long>::max()));
483  if(m_NumberOfBoundary == 0)
484  {
485  ok=0;
486  }
487  while( (m_NumberOfSeedsToAdded != 0) && ok)
488  {
489  count++;
490  m_VDGenerator->AddSeeds(m_NumberOfSeedsToAdded,m_SeedsToAdded.begin());
491  m_LastStepSeeds = m_NumberOfSeeds;
492  m_NumberOfSeeds += m_NumberOfSeedsToAdded;
493  this->RunSegmentOneStep();
494  // guess at a progress
495  this->UpdateProgress(static_cast<float>(count)
496  / static_cast<float>(NumericTraits<unsigned long>::max()));
497  }
498  }
499  else if(m_Steps == 1)
500  {
501  this->RunSegmentOneStep();
502  this->UpdateProgress(1.0);
503  }
504  else
505  {
506  this->RunSegmentOneStep();
507  if (m_Steps == 0)
508  {
509  this->UpdateProgress(1.0);
510  }
511  else
512  {
513  this->UpdateProgress(1.0/static_cast<float>(m_Steps));
514  }
515  if(m_NumberOfBoundary == 0)
516  {
517  ok=0;
518  }
519  int i = 1;
520  while((i<m_Steps) && ok)
521  {
522  m_VDGenerator->AddSeeds(m_NumberOfSeedsToAdded, m_SeedsToAdded.begin());
523  m_LastStepSeeds = m_NumberOfSeeds;
524  m_NumberOfSeeds += m_NumberOfSeedsToAdded;
525  this->RunSegmentOneStep();
526  i++;
527  this->UpdateProgress(static_cast<float>(i)/static_cast<float>(m_Steps));
528  }
529  }
530 }
531 
532 template <class TInputImage, class TOutputImage, class TBinaryPriorImage>
533 void
534 VoronoiSegmentationImageFilterBase <TInputImage,TOutputImage,TBinaryPriorImage>
535 ::GenerateData(void)
536 {
537  // Allocate the output
538  this->GetOutput()
539  ->SetBufferedRegion( this->GetOutput()->GetRequestedRegion() );
540  this->GetOutput()->Allocate();
541 
542  // This ivar should be necessary
543  m_Size = this->GetInput()->GetRequestedRegion().GetSize();
544 
545  VoronoiDiagram::PointType VDsize;
546  VDsize[0] = (VoronoiDiagram::CoordRepType)(m_Size[0]-0.1);
547  VDsize[1] = (VoronoiDiagram::CoordRepType)(m_Size[1]-0.1);
548  m_VDGenerator->SetBoundary(VDsize);
549  m_VDGenerator->SetRandomSeeds(m_NumberOfSeeds);
550 
551  this->RunSegment();
552  if(m_OutputBoundary)
553  {
554  this->MakeSegmentBoundary();
555  }
556  else
557  {
558  this->MakeSegmentObject();
559  }
560 
561 }
562 
563 template <class TInputImage, class TOutputImage, class TBinaryPriorImage>
564 void
565 VoronoiSegmentationImageFilterBase <TInputImage,TOutputImage,TBinaryPriorImage>
566 ::BeforeNextStep(void)
567 {
568  m_VDGenerator->AddSeeds(m_NumberOfSeedsToAdded, m_SeedsToAdded.begin());
569  m_LastStepSeeds = m_NumberOfSeeds;
570  m_NumberOfSeeds += m_NumberOfSeedsToAdded;
571 }
572 
573 
574 template <class TInputImage, class TOutputImage, class TBinaryPriorImage>
575 void
576 VoronoiSegmentationImageFilterBase <TInputImage,TOutputImage,TBinaryPriorImage>
577 ::MakeSegmentBoundary(void)
578 {
579 
580  RegionType region = this->GetInput()->GetRequestedRegion();
581  itk::ImageRegionIteratorWithIndex <OutputImageType> oit(this->GetOutput(), region);
582  while( !oit.IsAtEnd())
583  {
584  oit.Set(0);
585  ++oit;
586  }
587 
588  NeighborIdIterator nit;
589  NeighborIdIterator nitend;
590  for(int i=0;i<m_NumberOfSeeds;i++)
591  {
592  if(m_Label[i] == 2)
593  {
594  nitend = m_WorkingVD->NeighborIdsEnd(i);
595  for(nit=m_WorkingVD->NeighborIdsBegin(i); nit != nitend; ++nit)
596  {
597  if(((*nit)>i)&&(m_Label[*nit]==2))
598  {
599  drawLine(m_WorkingVD->GetSeed(i),m_WorkingVD->GetSeed(*nit));
600  i=i;
601  }
602  }
603  }
604  }
605 }
606 
607 template <class TInputImage, class TOutputImage, class TBinaryPriorImage>
608 void
609 VoronoiSegmentationImageFilterBase <TInputImage,TOutputImage,TBinaryPriorImage>
610 ::MakeSegmentObject(void)
611 {
612  RegionType region = this->GetInput()->GetRequestedRegion();
613  itk::ImageRegionIteratorWithIndex <OutputImageType> oit(this->GetOutput(), region);
614  while( !oit.IsAtEnd())
615  {
616  oit.Set(0);
617  ++oit;
618  }
619  PointIdIterator currPit;
620  PointIdIterator currPitEnd;
621  PointType currP;
622  PointTypeDeque VertList;
623  for(int i=0;i<m_NumberOfSeeds;i++)
624  {
625  if(m_Label[i] == 1)
626  {
627  CellAutoPointer currCell;
628  m_WorkingVD->GetCellId(i, currCell);
629  currPitEnd = currCell->PointIdsEnd();
630  VertList.clear();
631  for(currPit = currCell->PointIdsBegin(); currPit != currPitEnd; ++currPit)
632  {
633  m_WorkingVD->GetPoint((*currPit),&(currP));
634  VertList.push_back(currP);
635  }
636  FillPolygon(VertList);
637  }
638  }
639 }
640 
641 template <class TInputImage, class TOutputImage, class TBinaryPriorImage>
642 void
643 VoronoiSegmentationImageFilterBase <TInputImage,TOutputImage,TBinaryPriorImage>
644 ::FillPolygon(PointTypeDeque vertlist, OutputPixelType color)
645 {
646  IndexType idx;
647 
648  PointType currP;
649  PointType leftP;
650  PointType rightP;
651  currP = vertlist.front();
652  vertlist.pop_front();
653  leftP = vertlist.front();
654  while(currP[1] > leftP[1])
655  {
656  vertlist.push_back(currP);
657  currP=vertlist.front();
658  vertlist.pop_front();
659  leftP=vertlist.front();
660  }
661  rightP=vertlist.back();
662  while(currP[1] > rightP[1])
663  {
664  vertlist.push_front(currP);
665  currP=vertlist.back();
666  vertlist.pop_back();
667  rightP=vertlist.back();
668  }
669  leftP=vertlist.front();
670  PointTypeDeque tmpQ;
671  tmpQ.clear();
672  if(leftP[0]>rightP[0])
673  {
674  while(!(vertlist.empty()))
675  {
676  tmpQ.push_back(vertlist.front());
677  vertlist.pop_front();
678  }
679  while(!(tmpQ.empty()))
680  {
681  vertlist.push_front(tmpQ.front());
682  tmpQ.pop_front();
683  }
684  }
685  tmpQ.clear();
686  leftP=vertlist.front();
687  rightP=vertlist.back();
688 
689  double beginy=currP[1];
690  int intbeginy=(int)vcl_ceil(beginy);
691  idx[1]=intbeginy;
692  double leftendy=leftP[1];
693  double rightendy=rightP[1];
694  double beginx=currP[0];
695  double endx=currP[0];
696  double leftDx,rightDx;
697  double offset;
698  double leftheadx=beginx;
699  double rightheadx=endx;
700  double leftheady=beginy;
701  double rightheady=beginy;
702 
703  double endy;
704  bool RorL;
705  int i;
706  if(leftendy>rightendy)
707  {
708  RorL=1;
709  endy=rightendy;
710  }
711  else
712  {
713  RorL=0;
714  endy=leftendy;
715  }
716  if(leftP[1]==beginy)
717  {
718  leftDx = 1.0;
719  }
720  else
721  {
722  leftDx=(leftP[0]-beginx)/(leftP[1]-beginy);
723  }
724  if(rightP[1]==beginy)
725  {
726  rightDx = 1.0;
727  }
728  else
729  {
730  rightDx=(rightP[0]-endx)/(rightP[1]-beginy);
731  }
732  int intendy=(int)vcl_floor(endy);
733  if(intbeginy>intendy)
734  { //no scanline
735  if(RorL)
736  {
737  beginy=rightP[1];
738  }
739  else
740  {
741  beginy=leftP[1];
742  }
743  }
744  else if((intbeginy==intendy) && (intbeginy==0))
745  { //only one scanline at 0;
746  if(RorL)
747  {
748  endx=rightP[0];
749  }
750  else
751  {
752  beginx=leftP[0];
753  }
754  for(i=static_cast<int>(vcl_ceil(beginx));i<=static_cast<int>(vcl_floor(endx));i++)
755  {
756  idx[0]=i;
757  this->GetOutput()->SetPixel(idx,color);
758  }
759  idx[1] = idx[1] + 1;
760  }
761  else
762  { //normal case some scanlines
763  offset = (double)intbeginy-beginy;
764  endx += offset * rightDx;
765  beginx += offset * leftDx;
766  while(idx[1] <= intendy)
767  {
768  for(i = static_cast<int>(vcl_ceil(beginx)); i <= static_cast<int>(vcl_floor(endx)); i++)
769  {
770  idx[0] = i;
771  this->GetOutput()->SetPixel(idx,color);
772  }
773  endx += rightDx;
774  beginx += leftDx;
775  idx[1] = idx[1] + 1;
776  }
777  beginy=endy;
778  }
779 
780  int vsize = static_cast<int>( vertlist.size() );
781  while(vsize>2)
782  {
783  vsize--;
784  if(RorL)
785  {
786  vertlist.pop_back();
787  currP=rightP;
788  rightheadx=currP[0];
789  rightheady=currP[1];
790  endx=currP[0];
791  beginx=leftheadx+leftDx*(beginy-leftheady);
792  rightP=vertlist.back();
793  if (rightP[1] == currP[1])
794  {
795  rightDx = 1.0;
796  }
797  else
798  {
799  rightDx=(rightP[0]-currP[0])/(rightP[1]-currP[1]);
800  }
801  }
802  else
803  {
804  vertlist.pop_front();
805  currP=leftP;
806  leftheadx=currP[0];
807  leftheady=currP[1];
808  beginx=currP[0];
809  endx=rightheadx+rightDx*(beginy-rightheady);
810  leftP=vertlist.front();
811  if (leftP[1] == currP[1])
812  {
813  leftDx = 1.0;
814  }
815  else
816  {
817  leftDx=(leftP[0]-currP[0])/(leftP[1]-currP[1]);
818  }
819  }
820 
821  leftendy=leftP[1];
822  rightendy=rightP[1];
823  if(leftendy>rightendy)
824  {
825  RorL=1;
826  endy=rightendy;
827  }
828  else
829  {
830  RorL=0;
831  endy=leftendy;
832  }
833 
834  intendy=(int)vcl_floor(endy);
835  intbeginy=(int)vcl_ceil(beginy);
836 
837  if(intbeginy>intendy)
838  { //no scanline
839  if(RorL)
840  {
841  beginy=rightP[1];
842  }
843  else
844  {
845  beginy=leftP[1];
846  }
847  }
848  else
849  { //normal case some scanlines
850  offset = (double)intbeginy-beginy;
851  endx += offset * rightDx;
852  beginx += offset * leftDx;
853  while(idx[1] <= intendy)
854  {
855  for(i = static_cast<int>(vcl_ceil(beginx)); i <= static_cast<int>(vcl_floor(endx)); i++)
856  {
857  idx[0] = i;
858  this->GetOutput()->SetPixel(idx,color);
859  }
860  endx += rightDx;
861  beginx += leftDx;
862  idx[1] = idx[1] + 1;
863  }
864  beginy = idx[1];
865  }
866  }
867 
868 
869  if(RorL)
870  {
871  beginy=rightP[1];
872  endy=leftP[1];
873  }
874  else
875  {
876  beginy=leftP[1];
877  endy=rightP[1];
878  }
879  intbeginy=(int)vcl_ceil(beginy);
880  intendy=(int)vcl_floor(endy);
881  if(intbeginy<=intendy)
882  {
883  if(RorL)
884  {
885  if (rightP[1] == leftP[1])
886  {
887  rightDx = 1.0;
888  }
889  else
890  {
891  rightDx=(rightP[0]-leftP[0])/(rightP[1]-leftP[1]);
892  }
893  endx=rightP[0];
894  beginx=leftP[0]+leftDx*(rightP[1]-leftP[1]);
895  }
896  else
897  {
898  if (rightP[1] == leftP[1])
899  {
900  leftDx = 1.0;
901  }
902  else
903  {
904  leftDx=(rightP[0]-leftP[0])/(rightP[1]-leftP[1]);
905  }
906  beginx=leftP[0];
907  endx=rightP[0]+rightDx*(leftP[1]-rightP[1]);
908  }
909  offset = (double)intbeginy-beginy;
910  beginx += offset * leftDx;
911  endx += offset * rightDx;
912  while(idx[1] <= intendy)
913  {
914  for(i = static_cast<int>(vcl_ceil(beginx)); i <= static_cast<int>(vcl_floor(endx)); i++)
915  {
916  idx[0] = i;
917  this->GetOutput()->SetPixel(idx,color);
918  }
919  endx += rightDx;
920  beginx += leftDx;
921  idx[1] = idx[1] + 1;
922  }
923  }
924 
925 }
926 
927 template <class TInputImage, class TOutputImage, class TBinaryPriorImage>
928 void
929 VoronoiSegmentationImageFilterBase <TInputImage,TOutputImage,TBinaryPriorImage>
930 ::drawLine(PointType p1,PointType p2)
931 {
932  int x1=(int)(p1[0]+0.5);
933  int x2=(int)(p2[0]+0.5);
934  int y1=(int)(p1[1]+0.5);
935  int y2=(int)(p2[1]+0.5);
936  if(x1==static_cast<int>(m_Size[0])) x1--;
937  if(x2==static_cast<int>(m_Size[0])) x2--;
938  if(y1==static_cast<int>(m_Size[1])) y1--;
939  if(y2==static_cast<int>(m_Size[1])) y2--;
940 
941 
942  int dx=x1-x2;
943  int adx=(dx>0)?dx:-dx;
944  int dy=y1-y2;
945  int ady=(dy>0)?dy:-dy;
946  int save;
947  float curr;
948  IndexType idx;
949  if (adx > ady)
950  {
951  if(x1>x2)
952  {
953  save=x1; x1=x2; x2=save;
954  y1=y2;
955  }
956  curr=(float)y1;
957  if (dx == 0)
958  {
959  dx = 1;
960  }
961  float offset=(float)dy/dx;
962  for(int i=x1;i<=x2;i++)
963  {
964  idx[0]=i;
965  idx[1]=y1;
966  this->GetOutput()->SetPixel(idx,1);
967  curr += offset;
968  y1=(int)(curr+0.5);
969  }
970  }
971  else
972  {
973  if(y1>y2)
974  {
975  x1=x2;
976  save=y1; y1=y2; y2=save;
977  }
978  curr=(float)x1;
979  if (dy == 0)
980  {
981  dy = 1;
982  }
983  float offset=(float)dx/dy;
984  for(int i=y1;i<=y2;i++)
985  {
986  idx[0]=x1;
987  idx[1]=i;
988  this->GetOutput()->SetPixel(idx,1);
989  curr += offset;
990  x1=(int)(curr+0.5);
991  }
992  }
993 
994 }
995 
996 template <class TInputImage, class TOutputImage, class TBinaryPriorImage>
997 void
998 VoronoiSegmentationImageFilterBase <TInputImage,TOutputImage,TBinaryPriorImage>
999 ::DrawDiagram(VDImagePointer result,unsigned char incolor,
1000  unsigned char outcolor,unsigned char boundcolor)
1001 {
1002 
1003  RegionType region = this->GetInput()->GetRequestedRegion();
1004  itk::ImageRegionIteratorWithIndex <VDImage> vdit(result, region);
1005  while( !vdit.IsAtEnd())
1006  {
1007  vdit.Set(0);
1008  ++vdit;
1009  }
1010 
1011  EdgeIterator eit;
1012  EdgeIterator eitend = m_WorkingVD->EdgeEnd();
1013  PointType adds;
1014  Point<int,2> seeds;
1015  for(eit = m_WorkingVD->EdgeBegin();eit != eitend; ++eit)
1016  {
1017  seeds = m_WorkingVD->GetSeedsIDAroundEdge(&*eit);
1018  if((m_Label[seeds[0]]==2)||(m_Label[seeds[1]]==2))
1019  {
1020  drawVDline(result,eit->m_Left,eit->m_Right,boundcolor);
1021  }
1022  else if(m_Label[seeds[0]])
1023  {
1024  drawVDline(result,eit->m_Left,eit->m_Right,incolor);
1025  }
1026  else
1027  {
1028  drawVDline(result,eit->m_Left,eit->m_Right,outcolor);
1029  }
1030  }
1031 
1032 }
1033 
1034 template <class TInputImage, class TOutputImage, class TBinaryPriorImage>
1035 void
1036 VoronoiSegmentationImageFilterBase <TInputImage,TOutputImage,TBinaryPriorImage>
1037 ::drawVDline(VDImagePointer result,PointType p1,PointType p2,
1038  unsigned char color)
1039 {
1040  int x1=(int)(p1[0]+0.5);
1041  int x2=(int)(p2[0]+0.5);
1042  int y1=(int)(p1[1]+0.5);
1043  int y2=(int)(p2[1]+0.5);
1044  if(x1==(int)m_Size[0])
1045  {
1046  x1--;
1047  }
1048  if(x2==(int)m_Size[0])
1049  {
1050  x2--;
1051  }
1052  if(y1==(int)m_Size[1])
1053  {
1054  y1--;
1055  }
1056  if(y2==(int)m_Size[1])
1057  {
1058  y2--;
1059  }
1060  int dx=x1-x2;
1061  int adx=(dx>0)?dx:-dx;
1062  int dy=y1-y2;
1063  int ady=(dy>0)?dy:-dy;
1064  int save;
1065  float curr;
1066  IndexType idx;
1067  if (adx > ady)
1068  {
1069  if(x1>x2)
1070  {
1071  save=x1; x1=x2; x2=save;
1072  save=y1; y1=y2; y2=save;
1073  }
1074  curr=(float)y1;
1075  if (dx == 0)
1076  {
1077  dx = 1;
1078  }
1079  float offset=(float)dy/dx;
1080  for(int i=x1;i<=x2;i++)
1081  {
1082  idx[0]=i;
1083  idx[1]=y1;
1084  result->SetPixel(idx,color);
1085  curr += offset;
1086  y1=(int)(curr+0.5);
1087  }
1088  }
1089  else
1090  {
1091  if(y1>y2)
1092  {
1093  save=x1; x1=x2; x2=save;
1094  save=y1; y1=y2; y2=save;
1095  }
1096  curr=(float)x1;
1097  if (dy == 0)
1098  {
1099  dy = 1;
1100  }
1101  float offset=(float)dx/dy;
1102  for(int i=y1;i<=y2;i++)
1103  {
1104  idx[0]=x1;
1105  idx[1]=i;
1106  result->SetPixel(idx,color);
1107  curr += offset;
1108  x1=(int)(curr+0.5);
1109  }
1110  }
1111 }
1112 
1113 template <class TInputImage, class TOutputImage, class TBinaryPriorImage>
1114 void
1115 VoronoiSegmentationImageFilterBase <TInputImage,TOutputImage,TBinaryPriorImage>
1116 ::GenerateInputRequestedRegion()
1117 {
1118  // call the superclass's method
1119  Superclass::GenerateInputRequestedRegion();
1120 
1121  // set the input requested region to the LargestPossibleRegion
1122  {
1123  InputImagePointer input =
1124  const_cast< InputImageType * >( this->GetInput() );
1125  if (input)
1126  {
1127  input->SetRequestedRegion( this->GetInput()->GetLargestPossibleRegion() );
1128  }
1129  }
1130 }
1131 
1132 template <class TInputImage, class TOutputImage, class TBinaryPriorImage>
1133 void
1134 VoronoiSegmentationImageFilterBase <TInputImage,TOutputImage,TBinaryPriorImage>
1135 ::EnlargeOutputRequestedRegion(DataObject *output)
1136 {
1137  // call the superclass's method
1138  Superclass::EnlargeOutputRequestedRegion(output);
1139 
1140  // set the output requested region to the LargestPossibleRegion
1141  this->GetOutput()
1142  ->SetRequestedRegion( this->GetOutput()->GetLargestPossibleRegion() );
1143 }
1144 
1145 } //end namespace
1146 
1147 #endif
Generated at Sun May 19 2013 00:14:08 for Orfeo Toolbox with doxygen 1.8.3.1