GOFIGURE2  0.9.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
vtkFillImageWithPolyData.cxx
Go to the documentation of this file.
1 /*========================================================================
2  Copyright (c) INRIA - ASCLEPIOS Project (http://www-sop.inria.fr/asclepios).
3  All rights reserved.
4 
5  Redistribution and use in source and binary forms, with or without
6  modification, are permitted provided that the following conditions are met:
7 
8  * Redistributions of source code must retain the above copyright notice,
9  this list of conditions and the following disclaimer.
10 
11  * Redistributions in binary form must reproduce the above copyright notice,
12  this list of conditions and the following disclaimer in the documentation
13  and/or other materials provided with the distribution.
14 
15  * Neither the name of INRIA or ASCLEPIOS, nor the names of any contributors
16  may be used to endorse or promote products derived from this software
17  without specific prior written permission.
18 
19  * Modified source versions must be plainly marked as such, and must not be
20  misrepresented as being the original software.
21 
22  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND CONTRIBUTORS ``AS IS''
23  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
24  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
25  ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR
26  ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
27  DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
28  SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
29  CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
30  OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
31  OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32  =========================================================================*/
33 
34 /*=========================================================================
35  Modifications were made by the GoFigure Dev. Team.
36  while at Megason Lab, Systems biology, Harvard Medical school, 2009-11
37 
38  Copyright (c) 2009-11, President and Fellows of Harvard College.
39  All rights reserved.
40 
41  Redistribution and use in source and binary forms, with or without
42  modification, are permitted provided that the following conditions are met:
43 
44  Redistributions of source code must retain the above copyright notice,
45  this list of conditions and the following disclaimer.
46  Redistributions in binary form must reproduce the above copyright notice,
47  this list of conditions and the following disclaimer in the documentation
48  and/or other materials provided with the distribution.
49  Neither the name of the President and Fellows of Harvard College
50  nor the names of its contributors may be used to endorse or promote
51  products derived from this software without specific prior written
52  permission.
53 
54  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
55  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
56  THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
57  PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
58  BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
59  OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
60  OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
61  OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
62  WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
63  OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
64  ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
65 
66  =========================================================================*/
67 
69 
70 #include "vtkImageData.h"
71 #include "vtkInformation.h"
72 #include "vtkInformationVector.h"
73 #include "vtkObjectFactory.h"
74 #include "vtkStreamingDemandDrivenPipeline.h"
75 #include "vtkPoints.h"
76 
77 #include "vtkMath.h"
78 // extern int vtkrint(double a);
79 
80 vtkCxxRevisionMacro (vtkFillImageWithPolyData, "$Revision: 490 $");
82 
84 {
85  this->PolyData = 0;
86  this->InsidePixelValue = 1.0;
87  this->ExtractionDirection = 0;
88 }
89 
91 {
92  if ( this->PolyData ) { this->PolyData->Delete(); }
93 }
94 
95 //----------------------------------------------------------------------------
96 // The output extent is the intersection.
98  vtkInformation *vtkNotUsed(request),
99  vtkInformationVector **inputVector,
100  vtkInformationVector *outputVector)
101 {
102  (void)inputVector;
103 
104  // get the info objects
105  vtkInformation *outInfo = outputVector->GetInformationObject(0);
106 
107  // copy the polydata bounding box into bbox
108  if ( this->PolyData ) { this->PolyData->GetBounds (BBox); }
109 
110  vtkDataObject::SetPointDataActiveScalarInfo(outInfo, VTK_UNSIGNED_CHAR, 1);
111 
112  return 1;
113 }
114 
115 void vtkFillImageWithPolyData::PrintSelf(ostream & os, vtkIndent indent)
116 {
117  this->Superclass::PrintSelf(os, indent);
118  if ( this->PolyData )
119  {
120  os << indent << "PolyData: \n";
121  os << indent << *PolyData << endl;
122  }
123 }
124 
125 //----------------------------------------------------------------------------
126 // This templated function executes the filter for any type of data.
127 template< class T >
129  int ext[6],
130  vtkImageData *inData,
131  T *in1Ptr,
132  vtkImageData *outData,
133  unsigned char *outPtr,
134  int id, int slice, double bbox[6])
135 {
136  int idx0, idx1, idx2;
137  vtkIdType in1Inc0, in1Inc1, in1Inc2;
138  vtkIdType outInc0, outInc1, outInc2;
139  unsigned long count = 0;
140  unsigned long target;
141 
142 // numC = outData->GetNumberOfScalarComponents();
143 // pixSize = numC * static_cast< int >( sizeof( T ) );
144 
145  // Get information to march through data
146  inData->GetContinuousIncrements(ext, in1Inc0, in1Inc1, in1Inc2);
147  outData->GetContinuousIncrements(ext, outInc0, outInc1, outInc2);
148  int num0 = ext[1] - ext[0] + 1;
149  int num1 = ext[3] - ext[2] + 1;
150  int num2 = ext[5] - ext[4] + 1;
151 
152  double outVal = self->GetInsidePixelValue();
153 
154  vtkPoints *points = self->GetPolyData()->GetPoints();
155  vtkIdType numP = self->GetPolyData()->GetNumberOfPoints();
156 
157  double *spacing = inData->GetSpacing();
158  double *origin = inData->GetOrigin();
159 
160  target = (unsigned long)( num2 * num1 / 50.0 );
161  target++;
162 
163  // Loop through ouput pixels
164  for ( idx2 = ext[4]; idx2 < ext[4] + num2; ++idx2 )
165  {
166  for ( idx1 = ext[2]; !self->AbortExecute && idx1 < ext[2] + num1; ++idx1 )
167  {
168  if ( !id )
169  {
170  if ( !( count % target ) )
171  {
172  self->UpdateProgress( static_cast< double >( count )
173  / ( 50.0 * static_cast< double >( target ) ) );
174  }
175  count++;
176  }
177 
178  for ( idx0 = ext[0]; idx0 < ext[0] + num0; ++idx0 )
179  {
180  double val = 0.0;
181  if ( idx2 == slice )
182  {
183  double x = idx0 * spacing[0] + origin[0];
184  double y = idx1 * spacing[1] + origin[1];
185 
186  if ( x >= bbox[0] && x <= bbox[1] && y >= bbox[2] && y <= bbox[3] )
187  {
188  double angle = 0;
189 
190  for ( int i = 0; i < numP; i++ )
191  {
192  double p1[3], p2[3], dp1[2], dp2[2];
193  points->GetPoint (i, p1);
194  points->GetPoint ( ( i + 1 ) % numP, p2 );
195 
196  dp1[0] = p1[0] - x;
197  dp1[1] = p1[1] - y;
198 
199  dp2[0] = p2[0] - x;
200  dp2[1] = p2[1] - y;
201 
202  angle += self->Angle2D (dp1, dp2);
203  }
204  if ( fabs (angle) >= vtkMath::Pi() ) { val = outVal; }
205  }
206  }
207 
208  *outPtr = (unsigned char)( val );
209  ++outPtr;
210  in1Ptr += in1Inc0;
211  }
212  in1Ptr += in1Inc1;
213  outPtr += outInc1;
214  }
215  in1Ptr += in1Inc2;
216  outPtr += outInc2;
217  }
218 }
219 
220 //----------------------------------------------------------------------------
221 // This templated function executes the filter for any type of data.
222 template< class T >
224  vtkImageData *inData, T *in1Ptr,
225  vtkImageData *outData, unsigned char *outPtr,
226  int id, int slice, double bbox[6])
227 {
228  int idx0, idx1, idx2;
229  vtkIdType in1Inc0, in1Inc1, in1Inc2;
230  vtkIdType outInc0, outInc1, outInc2;
231  unsigned long count = 0;
232  unsigned long target;
233 
234 // numC = outData->GetNumberOfScalarComponents();
235 // int pixSize = numC * static_cast<int>(sizeof(T));
236 
237  // Get information to march through data
238  inData->GetContinuousIncrements(ext, in1Inc0, in1Inc1, in1Inc2);
239  outData->GetContinuousIncrements(ext, outInc0, outInc1, outInc2);
240  int num0 = ext[1] - ext[0] + 1;
241  int num1 = ext[3] - ext[2] + 1;
242  int num2 = ext[5] - ext[4] + 1;
243 
244  vtkPoints *points = self->GetPolyData()->GetPoints();
245  vtkIdType numP = self->GetPolyData()->GetNumberOfPoints();
246 
247  double *spacing = inData->GetSpacing();
248  double *origin = inData->GetOrigin();
249 
250  target = (unsigned long)( num2 * num1 / 50.0 );
251  target++;
252 
253  // Loop through ouput pixels
254  for ( idx2 = ext[4]; idx2 < ext[4] + num2; ++idx2 )
255  {
256  for ( idx1 = ext[2]; !self->AbortExecute && idx1 < ext[2] + num1; ++idx1 )
257  {
258  if ( !id )
259  {
260  if ( !( count % target ) )
261  {
262  self->UpdateProgress( static_cast< double >( count )
263  / ( 50.0 * static_cast< double >( target ) ) );
264  }
265  count++;
266  }
267 
268  for ( idx0 = ext[0]; idx0 < ext[0] + num0; ++idx0 )
269  {
270  double val = 0.0;
271 
272  if ( idx1 == slice )
273  {
274  double angle = 0;
275 
276  double x = idx0 * spacing[0] + origin[0];
277  double y = idx2 * spacing[2] + origin[2];
278 
279  if ( x >= bbox[0] && x <= bbox[1] && y >= bbox[4] && y <= bbox[5] )
280  {
281  for ( int i = 0; i < numP; i++ )
282  {
283  double p1[3], p2[3], dp1[2], dp2[2];
284  points->GetPoint (i, p1);
285  points->GetPoint ( ( i + 1 ) % numP, p2 );
286 
287  dp1[0] = p1[0] - x;
288  dp1[1] = p1[2] - y;
289 
290  dp2[0] = p2[0] - x;
291  dp2[1] = p2[2] - y;
292 
293  angle += self->Angle2D (dp1, dp2);
294  }
295 
296  if ( fabs (angle) >= vtkMath::Pi() ) { val = self->GetInsidePixelValue(); }
297  }
298  }
299 
300  *outPtr = (unsigned char)( val );
301  ++outPtr;
302  in1Ptr += in1Inc0;
303  }
304  in1Ptr += in1Inc1;
305  outPtr += outInc1;
306  }
307  in1Ptr += in1Inc2;
308  outPtr += outInc2;
309  }
310 }
311 
312 //----------------------------------------------------------------------------
313 // This templated function executes the filter for any type of data.
314 template< class T >
316  vtkImageData *inData, T *in1Ptr,
317  vtkImageData *outData, unsigned char *outPtr,
318  int id, int slice, double bbox[6])
319 {
320  int num0, num1, num2;
321  int idx0, idx1, idx2;
322  vtkIdType in1Inc0, in1Inc1, in1Inc2;
323  vtkIdType outInc0, outInc1, outInc2;
324  unsigned long count = 0;
325  unsigned long target;
326 
327 // numC = outData->GetNumberOfScalarComponents();
328 // pixSize = numC * static_cast< int >( sizeof( T ) );
329 
330  // Get information to march through data
331  inData->GetContinuousIncrements(ext, in1Inc0, in1Inc1, in1Inc2);
332  outData->GetContinuousIncrements(ext, outInc0, outInc1, outInc2);
333  num0 = ext[1] - ext[0] + 1;
334  num1 = ext[3] - ext[2] + 1;
335  num2 = ext[5] - ext[4] + 1;
336 
337  //get the polydata bounding box
338  //double bbox[6];
339  //self->GetPolyData()->GetBounds (bbox);
340 
341  vtkPoints *points = self->GetPolyData()->GetPoints();
342  vtkIdType numP = self->GetPolyData()->GetNumberOfPoints();
343 
344  double *spacing = inData->GetSpacing();
345  double *origin = inData->GetOrigin();
346 
347  target = (unsigned long)( num2 * num1 / 50.0 );
348  target++;
349 
350  // Loop through ouput pixels
351  for ( idx2 = ext[4]; idx2 < ext[4] + num2; ++idx2 )
352  {
353  for ( idx1 = ext[2]; !self->AbortExecute && idx1 < ext[2] + num1; ++idx1 )
354  {
355  if ( !id )
356  {
357  if ( !( count % target ) )
358  {
359  self->UpdateProgress( static_cast< double >( count )
360  / ( 50.0 * static_cast< double >( target ) ) );
361  }
362  count++;
363  }
364 
365  for ( idx0 = ext[0]; idx0 < ext[0] + num0; ++idx0 )
366  {
367  double val = 0.0;
368 
369  if ( idx0 == slice )
370  {
371  double angle = 0;
372 
373  double x = idx1 * spacing[1] + origin[1];
374  double y = idx2 * spacing[2] + origin[2];
375 
376  if ( x >= bbox[2] && x <= bbox[3] && y >= bbox[4] && y <= bbox[5] )
377  {
378  for ( int i = 0; i < numP; i++ )
379  {
380  double p1[3], p2[3], dp1[2], dp2[2];
381  points->GetPoint (i, p1);
382  points->GetPoint ( ( i + 1 ) % numP, p2 );
383 
384  dp1[0] = p1[1] - x;
385  dp1[1] = p1[2] - y;
386 
387  dp2[0] = p2[1] - x;
388  dp2[1] = p2[2] - y;
389 
390  angle += self->Angle2D (dp1, dp2);
391  }
392 
393  if ( fabs (angle) >= vtkMath::Pi() ) { val = self->GetInsidePixelValue(); }
394  }
395  }
396 
397  *outPtr = (unsigned char)( val );
398  ++outPtr;
399  in1Ptr += in1Inc0;
400  }
401  in1Ptr += in1Inc1;
402  outPtr += outInc1;
403  }
404  in1Ptr += in1Inc2;
405  outPtr += outInc2;
406  }
407 }
408 
409 //----------------------------------------------------------------------------
410 // This method is passed a input and output Datas, and executes the filter
411 // algorithm to fill the output from the inputs.
412 // It just executes a switch statement to call the correct function for
413 // the Datas data types.
415  vtkInformation *vtkNotUsed(request),
416  vtkInformationVector **vtkNotUsed(inputVector),
417  vtkInformationVector *vtkNotUsed(outputVector),
418  vtkImageData ***inData,
419  vtkImageData **outData,
420  int outExt[6], int id)
421 {
422  void *inPtr1;
423  void *outPtr;
424 
425  if ( !PolyData )
426  {
427  vtkErrorMacro("PolyData not set");
428  return;
429  }
430 
431  if ( PolyData->GetNumberOfPoints() < 3 )
432  {
433  vtkErrorMacro("PolyData must have more than 2 points.");
434  return;
435  }
436 
437  // check for the orientation and slice to draw the ROI
438  int slice = 0;
439  double *spacing = inData[0][0]->GetSpacing();
440  double *origin = inData[0][0]->GetOrigin();
441 
442  vtkPoints *points = PolyData->GetPoints();
443  double pt[3];
444  points->GetPoint (0, pt);
445 
446  if ( ExtractionDirection > 2 || ExtractionDirection < 0 )
447  {
448  vtkErrorMacro("Dont know extraction direction.");
449  }
450  else
451  {
452 // slice = int( vtkrint ( (pt[0]-origin[0])/spacing[0] ));
453  slice = static_cast< int >(
454  ( pt[ExtractionDirection] - origin[ExtractionDirection] )
455  / spacing[ExtractionDirection] );
456  }
457 
458  inPtr1 = inData[0][0]->GetScalarPointerForExtent(outExt);
459  outPtr = outData[0]->GetScalarPointerForExtent(outExt);
460 
461  if ( ExtractionDirection == 2 )
462  {
463  switch ( inData[0][0]->GetScalarType() )
464  {
465  vtkTemplateMacro( vtkFillImageWithPolyDataExecuteZ(this, outExt,
466  inData[0][0], (VTK_TT *)( inPtr1 ),
467  outData[0], (unsigned char *)( outPtr ), id, slice, BBox) );
468  default:
469  vtkErrorMacro(<< "Execute: Unknown ScalarType");
470  return;
471  }
472  }
473  else
474  {
475  if ( ExtractionDirection == 1 )
476  {
477  switch ( inData[0][0]->GetScalarType() )
478  {
479  vtkTemplateMacro( vtkFillImageWithPolyDataExecuteY(this, outExt,
480  inData[0][0], (VTK_TT *)( inPtr1 ),
481  outData[0], (unsigned char *)( outPtr ), id, slice, BBox) );
482  default:
483  vtkErrorMacro(<< "Execute: Unknown ScalarType");
484  return;
485  }
486  }
487  else
488  {
489  switch ( inData[0][0]->GetScalarType() )
490  {
491  vtkTemplateMacro( vtkFillImageWithPolyDataExecuteX(this, outExt,
492  inData[0][0], (VTK_TT *)( inPtr1 ),
493  outData[0], (unsigned char *)( outPtr ), id, slice, BBox) );
494  default:
495  vtkErrorMacro(<< "Execute: Unknown ScalarType");
496  return;
497  }
498  }
499  }
500 }
501 
502 double vtkFillImageWithPolyData::Angle2D(const double dp1[2], const double dp2[2])
503 {
504  double dt, t1, t2;
505 
506  t1 = atan2 (dp1[0], dp1[1]);
507  t2 = atan2 (dp2[0], dp2[1]);
508 
509  dt = t2 - t1;
510 
511  while ( dt > vtkMath::Pi() )
512  {
513  dt -= 2. * static_cast< double >( vtkMath::Pi() );
514  }
515  while ( dt < -vtkMath::Pi() )
516  {
517  dt += 2. * static_cast< double >( vtkMath::Pi() );
518  }
519 
520  return dt;
521 }
void vtkFillImageWithPolyDataExecuteX(vtkFillImageWithPolyData *self, int ext[6], vtkImageData *inData, T *in1Ptr, vtkImageData *outData, unsigned char *outPtr, int id, int slice, double bbox[6])
vtkCxxRevisionMacro(vtkFillImageWithPolyData,"$Revision: 490 $")
virtual void ThreadedRequestData(vtkInformation *vtkNotUsed(request), vtkInformationVector **vtkNotUsed(inputVector), vtkInformationVector *vtkNotUsed(outputVector), vtkImageData ***inData, vtkImageData **outData, int extent[6], int threadId)
void PrintSelf(ostream &os, vtkIndent indent)
void vtkFillImageWithPolyDataExecuteY(vtkFillImageWithPolyData *self, int ext[6], vtkImageData *inData, T *in1Ptr, vtkImageData *outData, unsigned char *outPtr, int id, int slice, double bbox[6])
void vtkFillImageWithPolyDataExecuteZ(vtkFillImageWithPolyData *self, int ext[6], vtkImageData *inData, T *in1Ptr, vtkImageData *outData, unsigned char *outPtr, int id, int slice, double bbox[6])
double Angle2D(const double dp1[2], const double dp2[2])
Compute the angle between the 2 input points.
vtkStandardNewMacro(vtkFillImageWithPolyData)
virtual int RequestInformation(vtkInformation *vtkNotUsed(request), vtkInformationVector **inputVector, vtkInformationVector *outputVector)