1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkROIStencilSource.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 "vtkROIStencilSource.h"
16
17 #include "vtkMath.h"
18 #include "vtkDataArray.h"
19 #include "vtkImageData.h"
20 #include "vtkImageStencilData.h"
21 #include "vtkInformation.h"
22 #include "vtkInformationVector.h"
23 #include "vtkObjectFactory.h"
24 #include "vtkStreamingDemandDrivenPipeline.h"
25
26 #include <math.h>
27
28 vtkStandardNewMacro(vtkROIStencilSource);
29
30 //----------------------------------------------------------------------------
vtkROIStencilSource()31 vtkROIStencilSource::vtkROIStencilSource()
32 {
33 this->SetNumberOfInputPorts(0);
34
35 this->Shape = vtkROIStencilSource::BOX;
36
37 this->Bounds[0] = 0.0;
38 this->Bounds[1] = 0.0;
39 this->Bounds[2] = 0.0;
40 this->Bounds[3] = 0.0;
41 this->Bounds[4] = 0.0;
42 this->Bounds[5] = 0.0;
43 }
44
45 //----------------------------------------------------------------------------
~vtkROIStencilSource()46 vtkROIStencilSource::~vtkROIStencilSource()
47 {
48 }
49
50 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)51 void vtkROIStencilSource::PrintSelf(ostream& os, vtkIndent indent)
52 {
53 this->Superclass::PrintSelf(os,indent);
54
55 os << indent << "Shape: " << this->GetShapeAsString() << "\n";
56 os << indent << "Bounds: " << this->Bounds[0] << " "
57 << this->Bounds[1] << " " << this->Bounds[2] << " "
58 << this->Bounds[3] << " " << this->Bounds[4] << " "
59 << this->Bounds[5] << "\n";
60 }
61
62 //----------------------------------------------------------------------------
GetShapeAsString()63 const char *vtkROIStencilSource::GetShapeAsString()
64 {
65 switch (this->Shape)
66 {
67 case vtkROIStencilSource::BOX:
68 return "Box";
69 case vtkROIStencilSource::ELLIPSOID:
70 return "Ellipsoid";
71 case vtkROIStencilSource::CYLINDERX:
72 return "CylinderX";
73 case vtkROIStencilSource::CYLINDERY:
74 return "CylinderY";
75 case vtkROIStencilSource::CYLINDERZ:
76 return "CylinderZ";
77 }
78 return "";
79 }
80
81 //----------------------------------------------------------------------------
82 // tolerance for stencil operations
83
84 #define VTK_STENCIL_TOL 7.62939453125e-06
85
86 //----------------------------------------------------------------------------
87 // Compute a reduced extent based on the Center and Size of the shape.
88 //
89 // Also returns the center and radius in voxel-index units.
vtkROIStencilSourceSubExtent(vtkROIStencilSource * self,const double origin[3],const double spacing[3],const int extent[6],int subextent[6],double icenter[3],double iradius[3])90 static void vtkROIStencilSourceSubExtent(
91 vtkROIStencilSource *self,
92 const double origin[3], const double spacing[3], const int extent[6],
93 int subextent[6], double icenter[3], double iradius[3])
94 {
95 double bounds[6];
96 self->GetBounds(bounds);
97
98 for (int i = 0; i < 3; i++)
99 {
100 icenter[i] = (0.5*(bounds[2*i] + bounds[2*i+1]) - origin[i])/spacing[i];
101 iradius[i] = 0.5*(bounds[2*i+1] - bounds[2*i])/spacing[i];
102
103 if (iradius[i] < 0)
104 {
105 iradius[i] = -iradius[i];
106 }
107 iradius[i] += VTK_STENCIL_TOL;
108
109 double emin = icenter[i] - iradius[i];
110 double emax = icenter[i] + iradius[i];
111
112 subextent[2*i] = extent[2*i];
113 subextent[2*i+1] = extent[2*i+1];
114
115 if (extent[2*i] < emin)
116 {
117 subextent[2*i] = VTK_INT_MAX;
118 if (extent[2*i+1] >= emin)
119 {
120 subextent[2*i] = vtkMath::Floor(emin) + 1;
121 }
122 }
123
124 if (extent[2*i+1] > emax)
125 {
126 subextent[2*i+1] = VTK_INT_MIN;
127 if (extent[2*i] <= emax)
128 {
129 subextent[2*i+1] = vtkMath::Floor(emax);
130 }
131 }
132 }
133 }
134
135 //----------------------------------------------------------------------------
vtkROIStencilSourceBox(vtkROIStencilSource * self,vtkImageStencilData * data,const int extent[6],const double origin[3],const double spacing[3])136 static int vtkROIStencilSourceBox(
137 vtkROIStencilSource *self, vtkImageStencilData *data,
138 const int extent[6], const double origin[3], const double spacing[3])
139 {
140 int subextent[6];
141 double icenter[3];
142 double iradius[3];
143
144 vtkROIStencilSourceSubExtent(self, origin, spacing, extent,
145 subextent, icenter, iradius);
146
147 // for keeping track of progress
148 unsigned long count = 0;
149 unsigned long target = static_cast<unsigned long>(
150 (subextent[5] - subextent[4] + 1)*
151 (subextent[3] - subextent[2] + 1)/50.0);
152 target++;
153
154 for (int idZ = subextent[4]; idZ <= subextent[5]; idZ++)
155 {
156 for (int idY = subextent[2]; idY <= subextent[3]; idY++)
157 {
158 if (count%target == 0)
159 {
160 self->UpdateProgress(count/(50.0*target));
161 }
162 count++;
163
164 int r1 = subextent[0];
165 int r2 = subextent[1];
166
167 if (r2 >= r1)
168 {
169 data->InsertNextExtent(r1, r2, idY, idZ);
170 }
171 } // for idY
172 } // for idZ
173
174 return 1;
175 }
176
177 //----------------------------------------------------------------------------
vtkROIStencilSourceEllipsoid(vtkROIStencilSource * self,vtkImageStencilData * data,const int extent[6],const double origin[3],const double spacing[3])178 static int vtkROIStencilSourceEllipsoid(
179 vtkROIStencilSource *self, vtkImageStencilData *data,
180 const int extent[6], const double origin[3], const double spacing[3])
181 {
182 int subextent[6];
183 double icenter[3];
184 double iradius[3];
185
186 vtkROIStencilSourceSubExtent(self, origin, spacing, extent,
187 subextent, icenter, iradius);
188
189 // for keeping track of progress
190 unsigned long count = 0;
191 unsigned long target = static_cast<unsigned long>(
192 (subextent[5] - subextent[4] + 1)*
193 (subextent[3] - subextent[2] + 1)/50.0);
194 target++;
195
196 for (int idZ = subextent[4]; idZ <= subextent[5]; idZ++)
197 {
198 double z = (idZ - icenter[2])/iradius[2];
199
200 for (int idY = subextent[2]; idY <= subextent[3]; idY++)
201 {
202 if (count%target == 0)
203 {
204 self->UpdateProgress(count/(50.0*target));
205 }
206 count++;
207
208 double y = (idY - icenter[1])/iradius[1];
209 double x2 = 1.0 - y*y - z*z;
210 if (x2 < 0)
211 {
212 continue;
213 }
214 double x = sqrt(x2);
215
216 int r1 = subextent[0];
217 int r2 = subextent[1];
218 double xmin = icenter[0] - x*iradius[0];
219 double xmax = icenter[0] + x*iradius[0];
220
221 if (r1 < xmin)
222 {
223 r1 = vtkMath::Floor(xmin) + 1;
224 }
225 if (r2 > xmax)
226 {
227 r2 = vtkMath::Floor(xmax);
228 }
229
230 if (r2 >= r1)
231 {
232 data->InsertNextExtent(r1, r2, idY, idZ);
233 }
234 } // for idY
235 } // for idZ
236
237 return 1;
238 }
239
240 //----------------------------------------------------------------------------
vtkROIStencilSourceCylinderX(vtkROIStencilSource * self,vtkImageStencilData * data,const int extent[6],const double origin[3],const double spacing[3])241 static int vtkROIStencilSourceCylinderX(
242 vtkROIStencilSource *self, vtkImageStencilData *data,
243 const int extent[6], const double origin[3], const double spacing[3])
244 {
245 int subextent[6];
246 double icenter[3];
247 double iradius[3];
248
249 vtkROIStencilSourceSubExtent(self, origin, spacing, extent,
250 subextent, icenter, iradius);
251
252 // for keeping track of progress
253 unsigned long count = 0;
254 unsigned long target = static_cast<unsigned long>(
255 (subextent[5] - subextent[4] + 1)*
256 (subextent[3] - subextent[2] + 1)/50.0);
257 target++;
258
259 for (int idZ = subextent[4]; idZ <= subextent[5]; idZ++)
260 {
261 double z = (idZ - icenter[2])/iradius[2];
262
263 for (int idY = subextent[2]; idY <= subextent[3]; idY++)
264 {
265 if (count%target == 0)
266 {
267 self->UpdateProgress(count/(50.0*target));
268 }
269 count++;
270
271 double y = (idY - icenter[1])/iradius[1];
272 if (y*y + z*z > 1.0)
273 {
274 continue;
275 }
276
277 int r1 = subextent[0];
278 int r2 = subextent[1];
279
280 if (r2 >= r1)
281 {
282 data->InsertNextExtent(r1, r2, idY, idZ);
283 }
284 } // for idY
285 } // for idZ
286
287 return 1;
288 }
289
290 //----------------------------------------------------------------------------
vtkROIStencilSourceCylinderY(vtkROIStencilSource * self,vtkImageStencilData * data,int extent[6],double origin[3],double spacing[3])291 static int vtkROIStencilSourceCylinderY(
292 vtkROIStencilSource *self, vtkImageStencilData *data,
293 int extent[6], double origin[3], double spacing[3])
294 {
295 int subextent[6];
296 double icenter[3];
297 double iradius[3];
298
299 vtkROIStencilSourceSubExtent(self, origin, spacing, extent,
300 subextent, icenter, iradius);
301
302 // for keeping track of progress
303 unsigned long count = 0;
304 unsigned long target = static_cast<unsigned long>(
305 (subextent[5] - subextent[4] + 1)*
306 (subextent[3] - subextent[2] + 1)/50.0);
307 target++;
308
309 for (int idZ = subextent[4]; idZ <= subextent[5]; idZ++)
310 {
311 double z = (idZ - icenter[2])/iradius[2];
312
313 for (int idY = subextent[2]; idY <= subextent[3]; idY++)
314 {
315 if (count%target == 0)
316 {
317 self->UpdateProgress(count/(50.0*target));
318 }
319 count++;
320
321 double x2 = 1.0 - z*z;
322 if (x2 < 0)
323 {
324 continue;
325 }
326 double x = sqrt(x2);
327
328 int r1 = subextent[0];
329 int r2 = subextent[1];
330 double xmin = icenter[0] - x*iradius[0];
331 double xmax = icenter[0] + x*iradius[0];
332
333 if (r1 < xmin)
334 {
335 r1 = vtkMath::Floor(xmin) + 1;
336 }
337 if (r2 > xmax)
338 {
339 r2 = vtkMath::Floor(xmax);
340 }
341
342 if (r2 >= r1)
343 {
344 data->InsertNextExtent(r1, r2, idY, idZ);
345 }
346 } // for idY
347 } // for idZ
348
349 return 1;
350 }
351
352 //----------------------------------------------------------------------------
vtkROIStencilSourceCylinderZ(vtkROIStencilSource * self,vtkImageStencilData * data,int extent[6],double origin[3],double spacing[3])353 static int vtkROIStencilSourceCylinderZ(
354 vtkROIStencilSource *self, vtkImageStencilData *data,
355 int extent[6], double origin[3], double spacing[3])
356 {
357 int subextent[6];
358 double icenter[3];
359 double iradius[3];
360
361 vtkROIStencilSourceSubExtent(self, origin, spacing, extent,
362 subextent, icenter, iradius);
363
364 // for keeping track of progress
365 unsigned long count = 0;
366 unsigned long target = static_cast<unsigned long>(
367 (subextent[5] - subextent[4] + 1)*
368 (subextent[3] - subextent[2] + 1)/50.0);
369 target++;
370
371 for (int idZ = subextent[4]; idZ <= subextent[5]; idZ++)
372 {
373 for (int idY = subextent[2]; idY <= subextent[3]; idY++)
374 {
375 if (count%target == 0)
376 {
377 self->UpdateProgress(count/(50.0*target));
378 }
379 count++;
380
381 double y = (idY - icenter[1])/iradius[1];
382 double x2 = 1.0 - y*y;
383 if (x2 < 0)
384 {
385 continue;
386 }
387 double x = sqrt(x2);
388
389 int r1 = subextent[0];
390 int r2 = subextent[1];
391 double xmin = icenter[0] - x*iradius[0];
392 double xmax = icenter[0] + x*iradius[0];
393
394 if (r1 < xmin)
395 {
396 r1 = vtkMath::Floor(xmin) + 1;
397 }
398 if (r2 > xmax)
399 {
400 r2 = vtkMath::Floor(xmax);
401 }
402
403 if (r2 >= r1)
404 {
405 data->InsertNextExtent(r1, r2, idY, idZ);
406 }
407 } // for idY
408 } // for idZ
409
410 return 1;
411 }
412
413
414 //----------------------------------------------------------------------------
RequestData(vtkInformation * request,vtkInformationVector ** inputVector,vtkInformationVector * outputVector)415 int vtkROIStencilSource::RequestData(
416 vtkInformation *request,
417 vtkInformationVector **inputVector,
418 vtkInformationVector *outputVector)
419 {
420 int extent[6];
421 double origin[3];
422 double spacing[3];
423 int result = 1;
424
425 this->Superclass::RequestData(request, inputVector, outputVector);
426
427 vtkInformation *outInfo = outputVector->GetInformationObject(0);
428 vtkImageStencilData *data = vtkImageStencilData::SafeDownCast(
429 outInfo->Get(vtkDataObject::DATA_OBJECT()));
430
431 outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), extent);
432 outInfo->Get(vtkDataObject::ORIGIN(), origin);
433 outInfo->Get(vtkDataObject::SPACING(), spacing);
434
435 switch (this->Shape)
436 {
437 case vtkROIStencilSource::BOX:
438 result = vtkROIStencilSourceBox(
439 this, data, extent, origin, spacing);
440 break;
441 case vtkROIStencilSource::ELLIPSOID:
442 result = vtkROIStencilSourceEllipsoid(
443 this, data, extent, origin, spacing);
444 break;
445 case vtkROIStencilSource::CYLINDERX:
446 result = vtkROIStencilSourceCylinderX(
447 this, data, extent, origin, spacing);
448 break;
449 case vtkROIStencilSource::CYLINDERY:
450 result = vtkROIStencilSourceCylinderY(
451 this, data, extent, origin, spacing);
452 break;
453 case vtkROIStencilSource::CYLINDERZ:
454 result = vtkROIStencilSourceCylinderZ(
455 this, data, extent, origin, spacing);
456 break;
457 }
458
459 return result;
460 }
461