GOFIGURE2  0.9.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
QGoSegmentationAlgo.cxx
Go to the documentation of this file.
1 /*=========================================================================
2  Authors: The GoFigure Dev. Team.
3  at Megason Lab, Systems biology, Harvard Medical school, 2009-11
4 
5  Copyright (c) 2009-11, President and Fellows of Harvard College.
6  All rights reserved.
7 
8  Redistribution and use in source and binary forms, with or without
9  modification, are permitted provided that the following conditions are met:
10 
11  Redistributions of source code must retain the above copyright notice,
12  this list of conditions and the following disclaimer.
13  Redistributions in binary form must reproduce the above copyright notice,
14  this list of conditions and the following disclaimer in the documentation
15  and/or other materials provided with the distribution.
16  Neither the name of the President and Fellows of Harvard College
17  nor the names of its contributors may be used to endorse or promote
18  products derived from this software without specific prior written
19  permission.
20 
21  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
23  THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
24  PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
25  BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
26  OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
27  OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
28  OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
29  WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
30  OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
31  ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32 
33 =========================================================================*/
34 #include "QGoSegmentationAlgo.h"
35 
36 // Extract ROI
37 #include "vtkExtractVOI.h"
38 
39 // extract polydata 2d
40 #include "vtkMarchingSquares.h"
41 //reorder contour
42 #include "vtkStripper.h"
43 #include "vtkCellArray.h"
44 // extract polydata 3d
45 #include "vtkContourFilter.h"
46 #include "vtkFeatureEdges.h"
47 #include "vtkFillHolesFilter.h"
48 #include "vtkPolyDataConnectivityFilter.h"
49 #include "vtkWindowedSincPolyDataFilter.h"
50 
51 #include "vtkMath.h"
52 
53 // Decimation 2d
54 #include "vtkPolylineDecimation.h"
55 // Decimation 3d
56 #include "vtkDecimatePro.h"
57 
58 
59 // test code
60 #include <assert.h>
61 
62 
64  : QObject( iParent ), m_Decimate(true),
65  m_NumberOfPoints(100), m_AlgoWidget(NULL)
66 {
67 }
68 //-------------------------------------------------------------------------
69 
70 //-------------------------------------------------------------------------
72 {
73 }
74 //-------------------------------------------------------------------------
75 
76 //-------------------------------------------------------------------------
78 {
79  return this->m_AlgoWidget;
80 }
81 //-------------------------------------------------------------------------
82 
83 //-------------------------------------------------------------------------
84 std::vector<vtkImageData*>
86  const std::vector<double>& iBounds,
87  const std::vector< vtkSmartPointer< vtkImageData > > & iImages)
88 {
89  assert( iBounds.size() == 6 );
90 
91  // vector to be returned
92  std::vector<vtkImageData*> listOfImages;
93 
94  //iterator on the images
95  std::vector< vtkSmartPointer< vtkImageData > >::const_iterator
96  it = iImages.begin();
97 
98  while( it != iImages.end())
99  {
100  listOfImages.push_back( VTKExtractROI(iBounds, *it) );
101  ++it;
102  }
103 
104  return listOfImages;
105 }
106 //-------------------------------------------------------------------------
107 
108 //-------------------------------------------------------------------------
109 vtkImageData*
111 VTKExtractROI(const std::vector<double>& iBounds,
112  const vtkSmartPointer< vtkImageData > & iImage)
113 {
114  // make sure there
115  assert( iBounds.size() == 6 );
116 
117  double org[3];
118  iImage->GetOrigin( org );
119 
120  double spacing[3];
121  iImage->GetSpacing( spacing );
122 
123  std::vector< int > bounds_idx(6);
124 
125  int k = 0;
126  for( int i = 0; i < 3; i++ )
127  {
128  bounds_idx[k] = vtkMath::Round( ( iBounds[k] - org[i] ) / spacing[i] );
129  ++k;
130  bounds_idx[k] = vtkMath::Round( ( iBounds[k] - org[i] ) / spacing[i] );
131  ++k;
132  }
133 
134  vtkSmartPointer<vtkExtractVOI> extractVOI =
135  vtkSmartPointer<vtkExtractVOI>::New();
136  extractVOI->SetInput( iImage );
137  extractVOI->SetVOI( bounds_idx[0], bounds_idx[1],
138  bounds_idx[2], bounds_idx[3],
139  bounds_idx[4], bounds_idx[5]);
140  extractVOI->Update();
141 
142  vtkImageData* temp_image = vtkImageData::New();
143  temp_image->DeepCopy( extractVOI->GetOutput() );
144 
145  return temp_image;
146 }
147 //-------------------------------------------------------------------------
148 
149 //-------------------------------------------------------------------------
150 std::vector<vtkPolyData*>
153  std::vector<vtkImageData*>& iInputImage,
154  const double & iThreshold)
155 {
156  // vector to be returned
157  std::vector<vtkPolyData*> listOfPolys;
158 
159  //iterator on the images
160  std::vector<vtkImageData*>::iterator it = iInputImage.begin();
161 
162  while( it != iInputImage.end())
163  {
164  vtkPolyData* extractedPolyData = ExtractPolyData(*it, iThreshold);
165  if( extractedPolyData )
166  {
167  listOfPolys.push_back( extractedPolyData );
168  }
169  ++it;
170  }
171 
172  return listOfPolys;
173 }
174 //-------------------------------------------------------------------------
175 
176 //-------------------------------------------------------------------------
177 vtkSmartPointer<vtkPolyData>
179 ExtractPolyData(vtkImageData *iInputImage, const double & iThreshold)
180 {
181  // make sure there is an image
182  assert( iInputImage );
183 
184  // test dimension of the input
185  int extent[6];
186  iInputImage->GetExtent( extent );
187 
188  int dimension = 3;
189 
190  for( int i = 0; i < 3; i++ )
191  {
192  if( extent[2*i] == extent[2*i+1] )
193  {
194  --dimension;
195  }
196  }
197 
198  switch ( dimension )
199  {
200  case 2 :
201  return ExtractContour( iInputImage, iThreshold);
202  case 3 :
203  return ExtractMesh( iInputImage, iThreshold);
204  default :
205  std::cout << "dimension unknown (Reconstruct polydata)" << std::endl;
206  }
207 
208  return NULL;
209 }
210 //-------------------------------------------------------------------------
211 
212 /*
213  * \todo Nicolas-do a better reconstruction..
214  */
215 //--------------------------------------------------------------------------
216 vtkSmartPointer<vtkPolyData>
218 ExtractContour( vtkSmartPointer<vtkImageData> iInputImage, const double & iThreshold)
219 {
220  // create iso-contours
221  vtkSmartPointer<vtkMarchingSquares> contours =
222  vtkSmartPointer<vtkMarchingSquares>::New();
223 
224  contours->SetInput(iInputImage);
225  contours->GenerateValues (1, iThreshold, iThreshold);
226  contours->Update();
227 
228  return ReorganizeContour( contours->GetOutput() );
229 }
230 //--------------------------------------------------------------------------
231 
232 //--------------------------------------------------------------------------
233 vtkSmartPointer<vtkPolyData>
235 ReorganizeContour(vtkSmartPointer<vtkPolyData> iInputImage)
236 {
237  /*
238  * \note Nicolas-to be enhanced
239  */
240  // Create reorganize contours
241  vtkStripper *stripper = vtkStripper::New();
242 
243  stripper->SetInput(iInputImage);
244  //Is it useful?? Which number is the best suited?
245  stripper->SetMaximumLength(999);
246  stripper->Update();
247 
248  // Reorder points
249  stripper->GetOutput()->GetLines()->InitTraversal();
250 
251  // npts = nb of points in the line
252  // *pts = pointer to each point
253 
254  vtkIdType *pts = NULL;
255  vtkIdType npts = 0;
256  stripper->GetOutput()->GetLines()->GetNextCell(npts, pts);
257  vtkPoints *points = vtkPoints::New();
258 
259  vtkCellArray *lines = vtkCellArray::New();
260  vtkIdType * lineIndices = new vtkIdType[static_cast< int >( npts + 1 )];
261 
262  for ( int k = 0; k < static_cast< int >( npts ); k++ )
263  {
264  points->InsertPoint( k, stripper->GetOutput()->GetPoints()->GetPoint(pts[k]) );
265  lineIndices[k] = k;
266  }
267 
268  lineIndices[static_cast< int >( npts )] = 0;
269  lines->InsertNextCell(npts + 1, lineIndices);
270  delete[] lineIndices;
271 
272  vtkSmartPointer<vtkPolyData> output = vtkSmartPointer<vtkPolyData>::New();
273  output->SetPoints(points);
274  output->SetLines(lines);
275 
276  lines->Delete();
277  points->Delete();
278  stripper->Delete();
279 
280  vtkPolyData* output2 = vtkPolyData::New();
281  output2->DeepCopy(output);
282 
283  return output2;
284 }
285 
286 //--------------------------------------------------------------------------
287 
288 //--------------------------------------------------------------------------
289 vtkSmartPointer<vtkPolyData>
291 ExtractMesh(vtkSmartPointer<vtkImageData> iInputImage, const double & iThreshold)
292 {
293  vtkSmartPointer< vtkContourFilter > contours = vtkSmartPointer< vtkContourFilter >::New();
294  contours->SetInput(iInputImage);
295  contours->SetComputeGradients(0);
296  contours->SetComputeNormals(0);
297  contours->SetComputeScalars(0);
298  contours->SetNumberOfContours(1);
299  contours->SetValue(0, iThreshold);
300  contours->Update();
301 
302  std::cout << "Number Of Points: " << iInputImage->GetNumberOfPoints()
303  << std::endl;
304 
305  vtkSmartPointer< vtkFeatureEdges > feature =
306  vtkSmartPointer< vtkFeatureEdges >::New();
307  feature->SetInputConnection( contours->GetOutputPort() );
308  feature->BoundaryEdgesOn();
309  feature->FeatureEdgesOff();
310  feature->NonManifoldEdgesOn();
311  feature->ManifoldEdgesOff();
312  feature->Update();
313 
314  vtkSmartPointer< vtkFillHolesFilter > fillFilter =
315  vtkSmartPointer< vtkFillHolesFilter >::New();
316 
317  vtkSmartPointer< vtkPolyDataConnectivityFilter > connectivityFilter =
318  vtkSmartPointer< vtkPolyDataConnectivityFilter >::New();
319  connectivityFilter->SetExtractionModeToLargestRegion();
320 
321 
322  if ( feature->GetOutput()->GetNumberOfCells() > 0 )
323  {
324  // fill holes if any!
325  fillFilter->SetInputConnection( contours->GetOutputPort() );
326  fillFilter->Update();
327 
328  connectivityFilter->SetInputConnection( fillFilter->GetOutputPort() );
329  }
330  else
331  {
332  connectivityFilter->SetInputConnection( contours->GetOutputPort() );
333  }
334 
335  // keep the largest region
336  connectivityFilter->Update();
337  /*
338  * \note Nicolas-Smoother not connected
339  */
340 /*
341  unsigned int smoothingIterations = 15;
342  double passBand = 0.001;
343  double featureAngle = 120.0;
344 
345  // smoothing
346  vtkSmartPointer< vtkWindowedSincPolyDataFilter > smoother =
347  vtkSmartPointer< vtkWindowedSincPolyDataFilter >::New();
348  smoother->SetInputConnection( connectivityFilter->GetOutputPort() );
349  smoother->SetNumberOfIterations( smoothingIterations );
350  smoother->BoundarySmoothingOff();
351  smoother->FeatureEdgeSmoothingOff();
352  smoother->SetFeatureAngle(featureAngle);
353  smoother->SetPassBand(passBand);
354  smoother->NonManifoldSmoothingOn();
355  smoother->NormalizeCoordinatesOn();
356  smoother->Update();
357 */
358  vtkPolyData* output = vtkPolyData::New();
359  output->DeepCopy( connectivityFilter->GetOutput() );
360 
361  return output;
362 }
363 //--------------------------------------------------------------------------
364 
365 //--------------------------------------------------------------------------
366 vtkSmartPointer<vtkPolyData>
369  vtkSmartPointer<vtkPolyData>& iPolyData,
370  const unsigned int& iNumberOfPoints)
371 {
372  // make sure there is a polydata
373  assert( iPolyData );
374 
375  /*
376  * \todo Nicolas-might be better solution to test dimension of a PolyData
377  */
378  // test dimension of the input
379  double* bounds = iPolyData->GetBounds();
380  int dimension = 3;
381 
382  for( int i = 0; i < 3; i++ )
383  {
384  if( bounds[2*i] == bounds[2*i+1] )
385  {
386  --dimension;
387  }
388  }
389 
390  switch ( dimension )
391  {
392  case 2 :
393  return DecimateContour( iPolyData, iNumberOfPoints);
394  case 3 :
395  return DecimateMesh( iPolyData, iNumberOfPoints);
396  default :
397  itkGenericExceptionMacro( << "dimension unknown (Reconstruct polydata)" );
398  }
399 
400  return NULL;
401 }
402 //--------------------------------------------------------------------------
403 
404 //--------------------------------------------------------------------------
405 vtkSmartPointer<vtkPolyData>
407 DecimateContour( vtkSmartPointer<vtkPolyData> iPolyData, unsigned int iNumberOfPoints)
408 {
409  // define target reduction
410  unsigned int numberOfPoints = iPolyData->GetNumberOfPoints();
411 
412  double target = 1 - (double)iNumberOfPoints/(double)numberOfPoints;
413 
414  vtkSmartPointer<vtkPolylineDecimation> decimator =
415  vtkSmartPointer<vtkPolylineDecimation>::New();
416  decimator->SetInput(iPolyData);
417  decimator->SetTargetReduction(target);
418  decimator->Update();
419 
420  return decimator->GetOutput();
421 }
422 //--------------------------------------------------------------------------
423 
424 //--------------------------------------------------------------------------
425 vtkSmartPointer<vtkPolyData>
427 DecimateMesh( vtkSmartPointer<vtkPolyData> iPolyData, unsigned int iNumberOfPoints)
428 {
429  // define target reduction
430  unsigned int numberOfPoints = iPolyData->GetNumberOfPoints();
431 
432  double target = 1 - (double)iNumberOfPoints/(double)numberOfPoints;
433 
434  vtkSmartPointer<vtkDecimatePro> decimator =
435  vtkSmartPointer<vtkDecimatePro>::New();
436  decimator->SetInput(iPolyData);
437  decimator->SetTargetReduction(target);
438  decimator->Update();
439 
440  /*
441  * \todo Nicolas-fill holes -smooth..?
442  */
443 
444  return decimator->GetOutput();
445 }
446 //--------------------------------------------------------------------------
447 
448 //--------------------------------------------------------------------------
449 void
451 SetDecimate(bool& iDecimate)
452 {
453  m_Decimate = iDecimate;
454 }
455 //--------------------------------------------------------------------------
456 
457 //--------------------------------------------------------------------------
458 bool
461 {
462  return m_Decimate;
463 }
464 //--------------------------------------------------------------------------
465 
466 //--------------------------------------------------------------------------
467 void
469 SetNumberOfPoints( const unsigned int& iNumberOfPoints)
470 {
471  m_NumberOfPoints = iNumberOfPoints;
472 }
473 //--------------------------------------------------------------------------
474 
475 //--------------------------------------------------------------------------
476 unsigned int
479 {
480  return m_NumberOfPoints;
481 }
482 //--------------------------------------------------------------------------
vtkSmartPointer< vtkPolyData > ReorganizeContour(vtkSmartPointer< vtkPolyData > iInputImage)
QGoAlgorithmWidget * m_AlgoWidget
vtkSmartPointer< vtkPolyData > DecimateContour(vtkSmartPointer< vtkPolyData > iPolyData, unsigned int iNumberOfPoints)
unsigned int GetNumberOfPoints() const
std::vector< vtkPolyData * > ExtractPolyData(std::vector< vtkImageData * > &iInputImage, const double &iThreshold)
std::vector< vtkImageData * > VTKExtractROI(const std::vector< double > &iBounds, const std::vector< vtkSmartPointer< vtkImageData > > &iImages)
return the vtkpolydata created by the algorithm
vtkSmartPointer< vtkPolyData > DecimatePolyData(vtkSmartPointer< vtkPolyData > &iPolyData, const unsigned int &iNumberOfPoints)
vtkSmartPointer< vtkPolyData > DecimateMesh(vtkSmartPointer< vtkPolyData > iPolyData, unsigned int iNumberOfPoints)
QGoAlgorithmWidget * GetAlgoWidget()
return the algowidget
void SetNumberOfPoints(const unsigned int &iNumberOfPoints)
QGoSegmentationAlgo(QWidget *iParent=0)
vtkSmartPointer< vtkPolyData > ExtractContour(vtkSmartPointer< vtkImageData > iInputImage, const double &iThreshold)
vtkSmartPointer< vtkPolyData > ExtractMesh(vtkSmartPointer< vtkImageData > iInputImage, const double &iThreshold)
void SetDecimate(bool &iDecimate)