1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkLassoStencilSource.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15 #include "vtkLassoStencilSource.h"
16
17 #include "vtkMath.h"
18 #include "vtkPoints.h"
19 #include "vtkCardinalSpline.h"
20 #include "vtkDataArray.h"
21 #include "vtkImageData.h"
22 #include "vtkImageStencilData.h"
23 #include "vtkInformation.h"
24 #include "vtkInformationVector.h"
25 #include "vtkObjectFactory.h"
26 #include "vtkStreamingDemandDrivenPipeline.h"
27 #include "vtkSmartPointer.h"
28
29 #include <cmath>
30 #include <map>
31
32 vtkStandardNewMacro(vtkLassoStencilSource);
33 vtkCxxSetObjectMacro(vtkLassoStencilSource, Points, vtkPoints);
34
35 //----------------------------------------------------------------------------
36 class vtkLSSPointMap : public std::map<int, vtkSmartPointer<vtkPoints> >
37 {
38 };
39
40 //----------------------------------------------------------------------------
vtkLassoStencilSource()41 vtkLassoStencilSource::vtkLassoStencilSource()
42 {
43 this->SetNumberOfInputPorts(0);
44
45 this->Shape = vtkLassoStencilSource::POLYGON;
46 this->SliceOrientation = 2;
47 this->Points = nullptr;
48 this->SplineX = vtkCardinalSpline::New();
49 this->SplineY = vtkCardinalSpline::New();
50
51 this->PointMap = new vtkLSSPointMap();
52 }
53
54 //----------------------------------------------------------------------------
~vtkLassoStencilSource()55 vtkLassoStencilSource::~vtkLassoStencilSource()
56 {
57 this->SetPoints(nullptr);
58 if (this->SplineX)
59 {
60 this->SplineX->Delete();
61 this->SplineX = nullptr;
62 }
63 if (this->SplineY)
64 {
65 this->SplineY->Delete();
66 this->SplineY = nullptr;
67 }
68 delete this->PointMap;
69 this->PointMap = nullptr;
70 }
71
72 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)73 void vtkLassoStencilSource::PrintSelf(ostream& os, vtkIndent indent)
74 {
75 this->Superclass::PrintSelf(os,indent);
76
77 os << indent << "Shape: " << this->GetShapeAsString() << "\n";
78 os << indent << "Points: " << this->Points << "\n";
79 os << indent << "SliceOrientation: " << this->GetSliceOrientation() << "\n";
80 os << indent << "SlicePoints: " << this->PointMap->size() << "\n";
81 }
82
83 //----------------------------------------------------------------------------
GetShapeAsString()84 const char *vtkLassoStencilSource::GetShapeAsString()
85 {
86 switch (this->Shape)
87 {
88 case vtkLassoStencilSource::POLYGON:
89 return "Polygon";
90 case vtkLassoStencilSource::SPLINE:
91 return "Spline";
92 }
93 return "";
94 }
95
96 //----------------------------------------------------------------------------
GetMTime()97 vtkMTimeType vtkLassoStencilSource::GetMTime()
98 {
99 vtkMTimeType mTime = this->vtkImageStencilSource::GetMTime();
100
101 if ( this->Points != nullptr )
102 {
103 vtkMTimeType t = this->Points->GetMTime();
104 if (t > mTime)
105 {
106 mTime = t;
107 }
108 }
109
110 if ( !this->PointMap->empty() )
111 {
112 vtkLSSPointMap::iterator iter = this->PointMap->begin();
113 while ( iter != this->PointMap->end() )
114 {
115 vtkMTimeType t = iter->second->GetMTime();
116 if (t > mTime)
117 {
118 mTime = t;
119 }
120 ++iter;
121 }
122 }
123
124 return mTime;
125 }
126
127 //----------------------------------------------------------------------------
SetSlicePoints(int i,vtkPoints * points)128 void vtkLassoStencilSource::SetSlicePoints(int i, vtkPoints *points)
129 {
130 vtkLSSPointMap::iterator iter = this->PointMap->find(i);
131 if (iter != this->PointMap->end())
132 {
133 if (iter->second == points)
134 {
135 return;
136 }
137 else if (points == nullptr)
138 {
139 this->PointMap->erase(iter);
140 }
141 else
142 {
143 iter->second = points;
144 }
145 }
146 else
147 {
148 if (points == nullptr)
149 {
150 return;
151 }
152 else
153 {
154 this->PointMap->insert(iter, vtkLSSPointMap::value_type(i, points));
155 }
156 }
157
158 this->Modified();
159 }
160
161 //----------------------------------------------------------------------------
RemoveAllSlicePoints()162 void vtkLassoStencilSource::RemoveAllSlicePoints()
163 {
164 this->PointMap->clear();
165 }
166
167 //----------------------------------------------------------------------------
GetSlicePoints(int i)168 vtkPoints *vtkLassoStencilSource::GetSlicePoints(int i)
169 {
170 vtkLSSPointMap::iterator iter = this->PointMap->find(i);
171 if (iter != this->PointMap->end())
172 {
173 return iter->second;
174 }
175 return nullptr;
176 }
177
178 //----------------------------------------------------------------------------
179 // tolerance for stencil operations
180
181 #define VTK_STENCIL_TOL 7.62939453125e-06
182
183 //----------------------------------------------------------------------------
184 // Compute a reduced extent based on the bounds of the shape.
vtkLassoStencilSourceSubExtent(vtkPoints * points,const double origin[3],const double spacing[3],const int extent[6],int subextent[6])185 static void vtkLassoStencilSourceSubExtent(
186 vtkPoints *points,
187 const double origin[3], const double spacing[3],
188 const int extent[6], int subextent[6])
189 {
190 double bounds[6];
191 points->GetBounds(bounds);
192
193 for (int i = 0; i < 3; i++)
194 {
195 double emin = (bounds[2*i] - origin[i])/spacing[i] - VTK_STENCIL_TOL;
196 double emax = (bounds[2*i+1] - origin[i])/spacing[i] + VTK_STENCIL_TOL;
197
198 subextent[2*i] = extent[2*i];
199 subextent[2*i+1] = extent[2*i+1];
200
201 if (extent[2*i] < emin)
202 {
203 subextent[2*i] = VTK_INT_MAX;
204 if (extent[2*i+1] >= emin)
205 {
206 subextent[2*i] = vtkMath::Floor(emin) + 1;
207 }
208 }
209
210 if (extent[2*i+1] > emax)
211 {
212 subextent[2*i+1] = VTK_INT_MIN;
213 if (extent[2*i] <= emax)
214 {
215 subextent[2*i+1] = vtkMath::Floor(emax);
216 }
217 }
218
219 }
220 }
221
222 //----------------------------------------------------------------------------
223 // Rasterize a polygon into the stencil
vtkLassoStencilSourcePolygon(vtkPoints * points,vtkImageStencilData * data,vtkImageStencilRaster * raster,const int extent[6],const double origin[3],const double spacing[3],int xj,int yj)224 static int vtkLassoStencilSourcePolygon(
225 vtkPoints *points, vtkImageStencilData *data, vtkImageStencilRaster *raster,
226 const int extent[6], const double origin[3], const double spacing[3],
227 int xj, int yj)
228 {
229 // get the bounds of the polygon
230 int subextent[6];
231 vtkLassoStencilSourceSubExtent(points, origin, spacing, extent, subextent);
232
233 // allocate the raster
234 raster->PrepareForNewData(&subextent[2*yj]);
235
236 // rasterize each line
237 vtkIdType n = points->GetNumberOfPoints();
238 double p[3];
239 double p0[2], p1[2], p2[2], p3[2];
240
241 points->GetPoint(n-1, p);
242 p0[0] = (p[xj] - origin[xj])/spacing[xj];
243 p0[1] = (p[yj] - origin[yj])/spacing[yj];
244
245 points->GetPoint(0, p);
246 p1[0] = (p[xj] - origin[xj])/spacing[xj];
247 p1[1] = (p[yj] - origin[yj])/spacing[yj];
248
249 double dx = p1[0] - p0[0];
250 double dy = p1[1] - p0[1];
251 if (dx*dx + dy*dy <= VTK_STENCIL_TOL*VTK_STENCIL_TOL)
252 {
253 n -= 1;
254 points->GetPoint(n-1, p);
255 p0[0] = (p[xj] - origin[xj])/spacing[xj];
256 p0[1] = (p[yj] - origin[yj])/spacing[yj];
257 }
258
259 points->GetPoint(1, p);
260 p2[0] = (p[xj] - origin[xj])/spacing[xj];
261 p2[1] = (p[yj] - origin[yj])/spacing[yj];
262
263 for (vtkIdType i = 0; i < n; i++)
264 {
265 points->GetPoint((i+2)%n, p);
266 p3[0] = (p[xj] - origin[xj])/spacing[xj];
267 p3[1] = (p[yj] - origin[yj])/spacing[yj];
268
269 raster->InsertLine(p1, p2);
270
271 p0[0] = p1[0]; p0[1] = p1[1];
272 p1[0] = p2[0]; p1[1] = p2[1];
273 p2[0] = p3[0]; p2[1] = p3[1];
274 }
275
276 raster->FillStencilData(data, extent, xj, yj);
277
278 return 1;
279 }
280
281
282 //----------------------------------------------------------------------------
283 // Generate the splines for the given set of points. The splines
284 // will be closed if the final point is equal to the first point.
285 // The parametric value for the resulting spline will be valid over
286 // the range [0, tmax] where the tmax value is returned by reference.
vtkLassoStencilSourceCreateSpline(vtkPoints * points,const double origin[3],const double spacing[3],int xj,int yj,vtkSpline * xspline,vtkSpline * yspline,double & tmax,double & dmax)287 static void vtkLassoStencilSourceCreateSpline(vtkPoints *points,
288 const double origin[3], const double spacing[3],
289 int xj, int yj, vtkSpline *xspline, vtkSpline *yspline,
290 double &tmax, double &dmax)
291 {
292 // initialize the spline
293 xspline->RemoveAllPoints();
294 yspline->RemoveAllPoints();
295 xspline->ClosedOff();
296 yspline->ClosedOff();
297
298 // get the number of points and the first/last point
299 vtkIdType n = points->GetNumberOfPoints();
300 double p[3];
301 double p0[2], p1[2];
302
303 points->GetPoint(n-1, p);
304 p0[0] = (p[xj] - origin[xj])/spacing[xj];
305 p0[1] = (p[yj] - origin[yj])/spacing[yj];
306
307 points->GetPoint(0, p);
308 p1[0] = (p[xj] - origin[xj])/spacing[xj];
309 p1[1] = (p[yj] - origin[yj])/spacing[yj];
310
311 // factor between real distance and parametric distance
312 double f = 1.0;
313 // the length of the implicit segment for closed loops
314 double lastd = 0;
315
316 // aspect ratio
317 double xf = 1.0;
318 double yf = 1.0;
319 if (spacing[xj] > spacing[yj])
320 {
321 xf = spacing[xj]/spacing[yj];
322 }
323 else
324 {
325 yf = spacing[yj]/spacing[xj];
326 }
327
328 // if first and last point are same, spline is closed
329 double dx = (p1[0] - p0[0])*xf;
330 double dy = (p1[1] - p0[1])*yf;
331 double d2 = dx*dx + dy*dy;
332 while (d2 <= VTK_STENCIL_TOL*VTK_STENCIL_TOL && n > 1)
333 {
334 n -= 1;
335 points->GetPoint(n-1, p);
336 p0[0] = (p[xj] - origin[xj])/spacing[xj];
337 p0[1] = (p[yj] - origin[yj])/spacing[yj];
338
339 xspline->ClosedOn();
340 yspline->ClosedOn();
341
342 // vtkSpline considers the parametric length of the implicit
343 // segment of closed loops to be unity, so set "f" so that
344 // multiplying the real length of that segment by "f" gives unity.
345 dx = (p1[0] - p0[0])*xf;
346 dy = (p1[1] - p0[1])*yf;
347 d2 = dx*dx + dy*dy;
348 lastd = sqrt(d2);
349 if (lastd > 0)
350 {
351 f = 1.0/lastd;
352 }
353 }
354
355 // Add all the points to the spline.
356 double d = 0.0;
357 for (vtkIdType i = 0; i < n; i++)
358 {
359 p0[0] = p1[0]; p0[1] = p1[1];
360
361 points->GetPoint(i, p);
362 p1[0] = (p[xj] - origin[xj])/spacing[xj];
363 p1[1] = (p[yj] - origin[yj])/spacing[yj];
364
365 dx = (p1[0] - p0[0])*xf;
366 dy = (p1[1] - p0[1])*yf;
367
368 d += sqrt(dx*dx + dy*dy);
369
370 double t = f*d;
371
372 xspline->AddPoint(t, p1[0]);
373 yspline->AddPoint(t, p1[1]);
374 }
375
376 // Do the spline precomputations
377 xspline->Compute();
378 yspline->Compute();
379
380 // The spline is valid over t = [0, tmax]
381 d += lastd;
382 tmax = f*d;
383 dmax = d;
384 }
385
386 //----------------------------------------------------------------------------
387 // Rasterize a spline contour into the stencil
vtkLassoStencilSourceSpline(vtkPoints * points,vtkImageStencilData * data,vtkImageStencilRaster * raster,const int extent[6],const double origin[3],const double spacing[3],int xj,int yj,vtkSpline * xspline,vtkSpline * yspline)388 static int vtkLassoStencilSourceSpline(
389 vtkPoints *points, vtkImageStencilData *data, vtkImageStencilRaster *raster,
390 const int extent[6], const double origin[3], const double spacing[3],
391 int xj, int yj, vtkSpline *xspline, vtkSpline *yspline)
392 {
393 // create the splines
394 double tmax, dmax;
395 vtkLassoStencilSourceCreateSpline(
396 points, origin, spacing, xj, yj, xspline, yspline, tmax, dmax);
397
398 if (dmax <= VTK_STENCIL_TOL)
399 {
400 return 1;
401 }
402
403 // get the bounds of the polygon as a first guess of the spline bounds
404 int subextent[6];
405 vtkLassoStencilSourceSubExtent(points, origin, spacing, extent, subextent);
406
407 // allocate the raster
408 raster->PrepareForNewData(&subextent[2*yj]);
409
410 // go around the spline
411 vtkIdType n = vtkMath::Floor(dmax)+1;
412 double delta = tmax/n;
413
414 double p1[2], p2[2], p3[2];
415
416 double t = 0;
417 p1[0] = xspline->Evaluate(t);
418 p1[1] = yspline->Evaluate(t);
419
420 t = delta;
421 p2[0] = xspline->Evaluate(t);
422 p2[1] = yspline->Evaluate(t);
423
424 for (vtkIdType i = 0; i < n; i++)
425 {
426 t += delta;
427 if (i == n-2)
428 {
429 t = 0;
430 }
431
432 p3[0] = xspline->Evaluate(t);
433 p3[1] = yspline->Evaluate(t);
434
435 raster->InsertLine(p1, p2);
436
437 p1[0] = p2[0]; p1[1] = p2[1];
438 p2[0] = p3[0]; p2[1] = p3[1];
439 }
440
441 raster->FillStencilData(data, extent, xj, yj);
442
443 return 1;
444 }
445
446 //----------------------------------------------------------------------------
vtkLassoStencilSourceExecute(vtkPoints * points,vtkImageStencilData * data,vtkImageStencilRaster * raster,int extent[6],double origin[3],double spacing[3],int shape,int xj,int yj,vtkSpline * xspline,vtkSpline * yspline)447 static int vtkLassoStencilSourceExecute(
448 vtkPoints *points, vtkImageStencilData *data, vtkImageStencilRaster *raster,
449 int extent[6], double origin[3], double spacing[3], int shape,
450 int xj, int yj, vtkSpline *xspline, vtkSpline *yspline)
451 {
452 int result = 1;
453
454 if (points == nullptr || points->GetNumberOfPoints() < 3)
455 {
456 return 1;
457 }
458
459 switch (shape)
460 {
461 case vtkLassoStencilSource::POLYGON:
462 result = vtkLassoStencilSourcePolygon(
463 points, data, raster, extent, origin, spacing, xj, yj);
464 break;
465 case vtkLassoStencilSource::SPLINE:
466 result = vtkLassoStencilSourceSpline(
467 points, data, raster, extent, origin, spacing, xj, yj,
468 xspline, yspline);
469 break;
470 }
471
472 return result;
473 }
474
475
476 //----------------------------------------------------------------------------
RequestData(vtkInformation * request,vtkInformationVector ** inputVector,vtkInformationVector * outputVector)477 int vtkLassoStencilSource::RequestData(
478 vtkInformation *request,
479 vtkInformationVector **inputVector,
480 vtkInformationVector *outputVector)
481 {
482 int extent[6];
483 double origin[3];
484 double spacing[3];
485 int result = 1;
486
487 this->Superclass::RequestData(request, inputVector, outputVector);
488
489 vtkInformation *outInfo = outputVector->GetInformationObject(0);
490 vtkImageStencilData *data = vtkImageStencilData::SafeDownCast(
491 outInfo->Get(vtkDataObject::DATA_OBJECT()));
492
493 outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), extent);
494 outInfo->Get(vtkDataObject::ORIGIN(), origin);
495 outInfo->Get(vtkDataObject::SPACING(), spacing);
496
497 int slabExtent[6];
498 slabExtent[0] = extent[0]; slabExtent[1] = extent[1];
499 slabExtent[2] = extent[2]; slabExtent[3] = extent[3];
500 slabExtent[4] = extent[4]; slabExtent[5] = extent[5];
501
502 int xj = 0;
503 int yj = 1;
504 int zj = 2;
505
506 if (this->SliceOrientation == 0)
507 {
508 xj = 1;
509 yj = 2;
510 zj = 0;
511 }
512 else if (this->SliceOrientation == 1)
513 {
514 xj = 0;
515 yj = 2;
516 zj = 1;
517 }
518
519 vtkImageStencilRaster raster(&extent[2*yj]);
520 raster.SetTolerance(VTK_STENCIL_TOL);
521
522 int zmin = extent[2*zj];
523 int zmax = extent[2*zj+1];
524
525 vtkLSSPointMap::iterator iter = this->PointMap->lower_bound(zmin);
526 vtkLSSPointMap::iterator maxiter = this->PointMap->upper_bound(zmax);
527
528 while (iter != maxiter && result != 0)
529 {
530 this->SetProgress((slabExtent[2*zj] - zmin)*1.0/(zmax - zmin + 1));
531
532 int i = iter->first;
533 vtkPoints *points = iter->second;
534
535 // fill in the slices with no SlicePoints
536 if (this->Points && i > slabExtent[2*zj])
537 {
538 slabExtent[2*zj+1] = i-1;
539
540 result = vtkLassoStencilSourceExecute(
541 this->Points, data, &raster, slabExtent, origin, spacing,
542 this->Shape, xj, yj, this->SplineX, this->SplineY);
543 }
544
545 // do the slice with its SlicePoints
546 if (result)
547 {
548 slabExtent[2*zj] = i;
549 slabExtent[2*zj+1] = i;
550
551 result = vtkLassoStencilSourceExecute(
552 points, data, &raster, slabExtent, origin, spacing,
553 this->Shape, xj, yj, this->SplineX, this->SplineY);
554
555 slabExtent[2*zj] = slabExtent[2*zj+1] + 1;
556 }
557
558 ++iter;
559 }
560
561 this->SetProgress((slabExtent[2*zj] - zmin)*1.0/(zmax - zmin + 1));
562
563 // fill in the rest
564 if (result && slabExtent[2*zj] <= zmax)
565 {
566 slabExtent[2*zj+1] = zmax;
567
568 result = vtkLassoStencilSourceExecute(
569 this->Points, data, &raster, slabExtent, origin, spacing,
570 this->Shape, xj, yj, this->SplineX, this->SplineY);
571
572 this->SetProgress(1.0);
573 }
574
575 return result;
576 }
577