1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkImageData.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
16 // Hide VTK_DEPRECATED_IN_9_0_0() warnings for this class.
17 #define VTK_DEPRECATION_LEVEL 0
18
19 #include "vtkImageData.h"
20
21 #include "vtkCellData.h"
22 #include "vtkDataArray.h"
23 #include "vtkGenericCell.h"
24 #include "vtkInformation.h"
25 #include "vtkInformationIntegerKey.h"
26 #include "vtkInformationVector.h"
27 #include "vtkLargeInteger.h"
28 #include "vtkLine.h"
29 #include "vtkMath.h"
30 #include "vtkMatrix3x3.h"
31 #include "vtkMatrix4x4.h"
32 #include "vtkObjectFactory.h"
33 #include "vtkPixel.h"
34 #include "vtkPointData.h"
35 #include "vtkPoints.h"
36 #include "vtkUnsignedCharArray.h"
37 #include "vtkVertex.h"
38 #include "vtkVoxel.h"
39
40 vtkStandardNewMacro(vtkImageData);
41 vtkStandardExtendedNewMacro(vtkImageData);
42
43 //------------------------------------------------------------------------------
vtkImageData()44 vtkImageData::vtkImageData()
45 {
46 int idx;
47
48 this->Vertex = nullptr;
49 this->Line = nullptr;
50 this->Pixel = nullptr;
51 this->Voxel = nullptr;
52
53 this->DataDescription = VTK_EMPTY;
54
55 for (idx = 0; idx < 3; ++idx)
56 {
57 this->Dimensions[idx] = 0;
58 this->Increments[idx] = 0;
59 this->Origin[idx] = 0.0;
60 this->Spacing[idx] = 1.0;
61 this->Point[idx] = 0.0;
62 }
63
64 this->DirectionMatrix = vtkMatrix3x3::New();
65 this->IndexToPhysicalMatrix = vtkMatrix4x4::New();
66 this->PhysicalToIndexMatrix = vtkMatrix4x4::New();
67 this->DirectionMatrix->Identity();
68 this->ComputeTransforms();
69
70 int extent[6] = { 0, -1, 0, -1, 0, -1 };
71 memcpy(this->Extent, extent, 6 * sizeof(int));
72
73 this->Information->Set(vtkDataObject::DATA_EXTENT_TYPE(), VTK_3D_EXTENT);
74 this->Information->Set(vtkDataObject::DATA_EXTENT(), this->Extent, 6);
75 }
76
77 //------------------------------------------------------------------------------
~vtkImageData()78 vtkImageData::~vtkImageData()
79 {
80 if (this->Vertex)
81 {
82 this->Vertex->Delete();
83 }
84 if (this->Line)
85 {
86 this->Line->Delete();
87 }
88 if (this->Pixel)
89 {
90 this->Pixel->Delete();
91 }
92 if (this->Voxel)
93 {
94 this->Voxel->Delete();
95 }
96 if (this->DirectionMatrix)
97 {
98 this->DirectionMatrix->Delete();
99 }
100 if (this->IndexToPhysicalMatrix)
101 {
102 this->IndexToPhysicalMatrix->Delete();
103 }
104 if (this->PhysicalToIndexMatrix)
105 {
106 this->PhysicalToIndexMatrix->Delete();
107 }
108 }
109
110 //------------------------------------------------------------------------------
111 // Copy the geometric and topological structure of an input structured points
112 // object.
CopyStructure(vtkDataSet * ds)113 void vtkImageData::CopyStructure(vtkDataSet* ds)
114 {
115 vtkImageData* sPts = static_cast<vtkImageData*>(ds);
116 this->Initialize();
117
118 int i;
119 for (i = 0; i < 3; i++)
120 {
121 this->Dimensions[i] = sPts->Dimensions[i];
122 this->Spacing[i] = sPts->Spacing[i];
123 this->Origin[i] = sPts->Origin[i];
124 }
125 this->DirectionMatrix->DeepCopy(sPts->GetDirectionMatrix());
126 this->ComputeTransforms();
127 this->SetExtent(sPts->GetExtent());
128 }
129
130 //------------------------------------------------------------------------------
Initialize()131 void vtkImageData::Initialize()
132 {
133 this->Superclass::Initialize();
134 if (this->Information)
135 {
136 this->SetDimensions(0, 0, 0);
137 }
138 }
139
140 //------------------------------------------------------------------------------
CopyInformationFromPipeline(vtkInformation * information)141 void vtkImageData::CopyInformationFromPipeline(vtkInformation* information)
142 {
143 // Let the superclass copy whatever it wants.
144 this->Superclass::CopyInformationFromPipeline(information);
145
146 // Copy origin and spacing from pipeline information to the internal
147 // copies.
148 if (information->Has(SPACING()))
149 {
150 this->SetSpacing(information->Get(SPACING()));
151 }
152 if (information->Has(ORIGIN()))
153 {
154 this->SetOrigin(information->Get(ORIGIN()));
155 }
156 if (information->Has(DIRECTION()))
157 {
158 this->SetDirectionMatrix(information->Get(DIRECTION()));
159 }
160 }
161
162 //------------------------------------------------------------------------------
CopyInformationToPipeline(vtkInformation * info)163 void vtkImageData::CopyInformationToPipeline(vtkInformation* info)
164 {
165 // Let the superclass copy information to the pipeline.
166 this->Superclass::CopyInformationToPipeline(info);
167
168 // Copy the spacing, origin, direction, and scalar info
169 info->Set(vtkDataObject::SPACING(), this->Spacing, 3);
170 info->Set(vtkDataObject::ORIGIN(), this->Origin, 3);
171 info->Set(vtkDataObject::DIRECTION(), this->DirectionMatrix->GetData(), 9);
172 vtkDataObject::SetPointDataActiveScalarInfo(
173 info, this->GetScalarType(), this->GetNumberOfScalarComponents());
174 }
175
176 //------------------------------------------------------------------------------
177 // Graphics filters reallocate every execute. Image filters try to reuse
178 // the scalars.
PrepareForNewData()179 void vtkImageData::PrepareForNewData()
180 {
181 // free everything but the scalars
182 vtkDataArray* scalars = this->GetPointData()->GetScalars();
183 if (scalars)
184 {
185 scalars->Register(this);
186 }
187 this->Initialize();
188 if (scalars)
189 {
190 this->GetPointData()->SetScalars(scalars);
191 scalars->UnRegister(this);
192 }
193 }
194
195 //------------------------------------------------------------------------------
GetCellDims(int cellDims[3])196 void vtkImageData::GetCellDims(int cellDims[3])
197 {
198 for (int i = 0; i < 3; ++i)
199 {
200 cellDims[i] = ((this->Dimensions[i] - 1) < 1) ? 1 : this->Dimensions[i] - 1;
201 }
202 }
203
204 namespace
205 {
206 class CellVisibility
207 {
208 public:
CellVisibility(vtkImageData * input)209 CellVisibility(vtkImageData* input)
210 : Input(input)
211 {
212 }
operator ()(const vtkIdType id)213 bool operator()(const vtkIdType id) { return !Input->IsCellVisible(id); }
214
215 private:
216 vtkImageData* Input;
217 };
218 } // anonymous namespace
219
220 //------------------------------------------------------------------------------
GetCellNeighbors(vtkIdType cellId,vtkIdList * ptIds,vtkIdList * cellIds)221 void vtkImageData::GetCellNeighbors(vtkIdType cellId, vtkIdList* ptIds, vtkIdList* cellIds)
222 {
223 int numPtIds = ptIds->GetNumberOfIds();
224
225 // Use special methods for speed
226 switch (numPtIds)
227 {
228 case 0:
229 cellIds->Reset();
230 return;
231
232 case 1:
233 case 2:
234 case 4: // vertex, edge, face neighbors
235 vtkStructuredData::GetCellNeighbors(cellId, ptIds, cellIds, this->GetDimensions());
236 break;
237
238 default:
239 this->Superclass::GetCellNeighbors(cellId, ptIds, cellIds);
240 }
241
242 // If blanking, remove blanked cells.
243 if (this->GetPointGhostArray() || this->GetCellGhostArray())
244 {
245 vtkIdType* pCellIds = cellIds->GetPointer(0);
246 vtkIdType* end =
247 std::remove_if(pCellIds, pCellIds + cellIds->GetNumberOfIds(), CellVisibility(this));
248 cellIds->Resize(std::distance(pCellIds, end));
249 }
250 }
251
252 //------------------------------------------------------------------------------
GetCellNeighbors(vtkIdType cellId,vtkIdList * ptIds,vtkIdList * cellIds,int * seedLoc)253 void vtkImageData::GetCellNeighbors(
254 vtkIdType cellId, vtkIdList* ptIds, vtkIdList* cellIds, int* seedLoc)
255 {
256 int numPtIds = ptIds->GetNumberOfIds();
257
258 // Use special methods for speed
259 switch (numPtIds)
260 {
261 case 0:
262 cellIds->Reset();
263 return;
264
265 case 1:
266 case 2:
267 case 4: // vertex, edge, face neighbors
268 vtkStructuredData::GetCellNeighbors(cellId, ptIds, cellIds, this->GetDimensions(), seedLoc);
269 break;
270
271 default:
272 this->Superclass::GetCellNeighbors(cellId, ptIds, cellIds);
273 }
274
275 // If blanking, remove blanked cells.
276 if (this->GetPointGhostArray() || this->GetCellGhostArray())
277 {
278 vtkIdType* pCellIds = cellIds->GetPointer(0);
279 vtkIdType* end =
280 std::remove_if(pCellIds, pCellIds + cellIds->GetNumberOfIds(), CellVisibility(this));
281 cellIds->Resize(std::distance(pCellIds, end));
282 }
283 }
284
285 //------------------------------------------------------------------------------
IsPointVisible(vtkIdType pointId)286 unsigned char vtkImageData::IsPointVisible(vtkIdType pointId)
287 {
288 return vtkStructuredData::IsPointVisible(pointId, this->GetPointGhostArray());
289 }
290
291 //------------------------------------------------------------------------------
292 // Return non-zero if the specified cell is visible (i.e., not blanked)
IsCellVisible(vtkIdType cellId)293 unsigned char vtkImageData::IsCellVisible(vtkIdType cellId)
294 {
295 return vtkStructuredData::IsCellVisible(cellId, this->Dimensions, this->DataDescription,
296 this->GetCellGhostArray(), this->GetPointGhostArray());
297 }
298
299 //------------------------------------------------------------------------------
300 template <class T>
vtkImageDataGetTypeSize(T *)301 unsigned long vtkImageDataGetTypeSize(T*)
302 {
303 return sizeof(T);
304 }
305
306 //------------------------------------------------------------------------------
GetCellTemplateForDataDescription()307 vtkCell* vtkImageData::GetCellTemplateForDataDescription()
308 {
309 vtkCell* cell = nullptr;
310 switch (this->DataDescription)
311 {
312 case VTK_EMPTY:
313 break;
314
315 case VTK_SINGLE_POINT:
316 cell = this->Vertex;
317 break;
318
319 case VTK_X_LINE:
320 case VTK_Y_LINE:
321 case VTK_Z_LINE:
322 cell = this->Line;
323 break;
324
325 case VTK_XY_PLANE:
326 case VTK_YZ_PLANE:
327 case VTK_XZ_PLANE:
328 cell = this->Pixel;
329 break;
330
331 case VTK_XYZ_GRID:
332 cell = this->Voxel;
333 break;
334
335 default:
336 vtkErrorMacro("Invalid DataDescription.");
337 break;
338 }
339 return cell;
340 }
341
342 //------------------------------------------------------------------------------
GetCellTemplateForDataDescription(vtkGenericCell * cell)343 bool vtkImageData::GetCellTemplateForDataDescription(vtkGenericCell* cell)
344 {
345 switch (this->DataDescription)
346 {
347 case VTK_EMPTY:
348 cell->SetCellTypeToEmptyCell();
349 break;
350
351 case VTK_SINGLE_POINT:
352 cell->SetCellTypeToVertex();
353 break;
354
355 case VTK_X_LINE:
356 case VTK_Y_LINE:
357 case VTK_Z_LINE:
358 cell->SetCellTypeToLine();
359 break;
360
361 case VTK_XY_PLANE:
362 case VTK_YZ_PLANE:
363 case VTK_XZ_PLANE:
364 cell->SetCellTypeToPixel();
365 break;
366
367 case VTK_XYZ_GRID:
368 cell->SetCellTypeToVoxel();
369 break;
370
371 default:
372 vtkErrorMacro("Invalid DataDescription.");
373 return false;
374 }
375 return true;
376 }
377
378 //------------------------------------------------------------------------------
GetIJKMinForCellId(vtkIdType cellId,int ijkMin[3])379 bool vtkImageData::GetIJKMinForCellId(vtkIdType cellId, int ijkMin[3])
380 {
381 vtkIdType dims[3];
382 this->GetDimensions(dims);
383
384 ijkMin[0] = ijkMin[1] = ijkMin[2] = 0;
385
386 if (dims[0] == 0 || dims[1] == 0 || dims[2] == 0)
387 {
388 vtkErrorMacro("Requesting a cell from an empty image.");
389 return false;
390 }
391
392 switch (this->DataDescription)
393 {
394 case VTK_EMPTY:
395 return false;
396
397 case VTK_SINGLE_POINT:
398 // cellId can only be = 0
399 break;
400
401 case VTK_X_LINE:
402 ijkMin[0] = cellId;
403 break;
404
405 case VTK_Y_LINE:
406 ijkMin[1] = cellId;
407 break;
408
409 case VTK_Z_LINE:
410 ijkMin[2] = cellId;
411 break;
412
413 case VTK_XY_PLANE:
414 ijkMin[0] = cellId % (dims[0] - 1);
415 ijkMin[1] = cellId / (dims[0] - 1);
416 break;
417
418 case VTK_YZ_PLANE:
419 ijkMin[1] = cellId % (dims[1] - 1);
420 ijkMin[2] = cellId / (dims[1] - 1);
421 break;
422
423 case VTK_XZ_PLANE:
424 ijkMin[0] = cellId % (dims[0] - 1);
425 ijkMin[2] = cellId / (dims[0] - 1);
426 break;
427
428 case VTK_XYZ_GRID:
429 ijkMin[0] = cellId % (dims[0] - 1);
430 ijkMin[1] = (cellId / (dims[0] - 1)) % (dims[1] - 1);
431 ijkMin[2] = cellId / ((dims[0] - 1) * (dims[1] - 1));
432 break;
433
434 default:
435 vtkErrorMacro("Invalid DataDescription.");
436 return false;
437 }
438 return true;
439 }
440
441 //------------------------------------------------------------------------------
GetIJKMaxForIJKMin(int ijkMin[3],int ijkMax[3])442 bool vtkImageData::GetIJKMaxForIJKMin(int ijkMin[3], int ijkMax[3])
443 {
444 vtkIdType dims[3];
445 this->GetDimensions(dims);
446
447 ijkMax[0] = ijkMax[1] = ijkMax[2] = 0;
448
449 if (dims[0] == 0 || dims[1] == 0 || dims[2] == 0)
450 {
451 vtkErrorMacro("Requesting a cell from an empty image.");
452 return false;
453 }
454
455 switch (this->DataDescription)
456 {
457 case VTK_EMPTY:
458 return false;
459
460 case VTK_SINGLE_POINT:
461 // cellId can only be = 0
462 break;
463
464 case VTK_X_LINE:
465 ijkMax[0] = ijkMin[0] + 1;
466 break;
467
468 case VTK_Y_LINE:
469 ijkMax[1] = ijkMin[1] + 1;
470 break;
471
472 case VTK_Z_LINE:
473 ijkMax[2] = ijkMin[2] + 1;
474 break;
475
476 case VTK_XY_PLANE:
477 ijkMax[0] = ijkMin[0] + 1;
478 ijkMax[1] = ijkMin[1] + 1;
479 break;
480
481 case VTK_YZ_PLANE:
482 ijkMax[1] = ijkMin[1] + 1;
483 ijkMax[2] = ijkMin[2] + 1;
484 break;
485
486 case VTK_XZ_PLANE:
487 ijkMax[0] = ijkMin[0] + 1;
488 ijkMax[2] = ijkMin[2] + 1;
489 break;
490
491 case VTK_XYZ_GRID:
492 ijkMax[0] = ijkMin[0] + 1;
493 ijkMax[1] = ijkMin[1] + 1;
494 ijkMax[2] = ijkMin[2] + 1;
495 break;
496
497 default:
498 vtkErrorMacro("Invalid DataDescription.");
499 return false;
500 }
501 return true;
502 }
503
504 //------------------------------------------------------------------------------
AddPointsToCellTemplate(vtkCell * cell,int ijkMin[3],int ijkMax[3])505 void vtkImageData::AddPointsToCellTemplate(vtkCell* cell, int ijkMin[3], int ijkMax[3])
506 {
507 int loc[3], i, j, k;
508 vtkIdType idx, npts;
509 double xyz[3];
510 const int* extent = this->Extent;
511
512 vtkIdType dims[3];
513 this->GetDimensions(dims);
514 vtkIdType d01 = dims[0] * dims[1];
515
516 // Extract point coordinates and point ids
517 // Ids are relative to extent min.
518 npts = 0;
519 for (loc[2] = ijkMin[2]; loc[2] <= ijkMax[2]; loc[2]++)
520 {
521 k = loc[2] + extent[4];
522 for (loc[1] = ijkMin[1]; loc[1] <= ijkMax[1]; loc[1]++)
523 {
524 j = loc[1] + extent[2];
525 for (loc[0] = ijkMin[0]; loc[0] <= ijkMax[0]; loc[0]++)
526 {
527 i = loc[0] + extent[0];
528 this->TransformIndexToPhysicalPoint(i, j, k, xyz);
529
530 idx = loc[0] + loc[1] * dims[0] + loc[2] * d01;
531 cell->PointIds->SetId(npts, idx);
532 cell->Points->SetPoint(npts++, xyz);
533 }
534 }
535 }
536 }
537
538 //------------------------------------------------------------------------------
GetCell(vtkIdType cellId)539 vtkCell* vtkImageData::GetCell(vtkIdType cellId)
540 {
541 int ijkMin[3];
542 if (!this->GetIJKMinForCellId(cellId, ijkMin))
543 {
544 return nullptr;
545 }
546
547 // Need to use vtImageData:: to avoid calling child classes implementation
548 return this->vtkImageData::GetCell(ijkMin[0], ijkMin[1], ijkMin[2]);
549 }
550
551 //------------------------------------------------------------------------------
GetCell(int iMin,int jMin,int kMin)552 vtkCell* vtkImageData::GetCell(int iMin, int jMin, int kMin)
553 {
554 vtkCell* cell = this->GetCellTemplateForDataDescription();
555 if (cell == nullptr)
556 {
557 return nullptr;
558 }
559
560 int ijkMin[3] = { iMin, jMin, kMin };
561 int ijkMax[3];
562 if (!this->GetIJKMaxForIJKMin(ijkMin, ijkMax))
563 {
564 return nullptr;
565 }
566
567 this->AddPointsToCellTemplate(cell, ijkMin, ijkMax);
568 return cell;
569 }
570
571 //------------------------------------------------------------------------------
GetCell(vtkIdType cellId,vtkGenericCell * cell)572 void vtkImageData::GetCell(vtkIdType cellId, vtkGenericCell* cell)
573 {
574 if (!this->GetCellTemplateForDataDescription(cell))
575 {
576 cell->SetCellTypeToEmptyCell();
577 return;
578 }
579
580 int ijkMin[3];
581 if (!this->GetIJKMinForCellId(cellId, ijkMin))
582 {
583 cell->SetCellTypeToEmptyCell();
584 return;
585 }
586
587 int ijkMax[3];
588 if (!this->GetIJKMaxForIJKMin(ijkMin, ijkMax))
589 {
590 cell->SetCellTypeToEmptyCell();
591 return;
592 }
593
594 this->AddPointsToCellTemplate(cell, ijkMin, ijkMax);
595 }
596
597 //------------------------------------------------------------------------------
598 // Fast implementation of GetCellBounds(). Bounds are calculated without
599 // constructing a cell.
GetCellBounds(vtkIdType cellId,double bounds[6])600 void vtkImageData::GetCellBounds(vtkIdType cellId, double bounds[6])
601 {
602 int ijkMin[3];
603 if (!this->GetIJKMinForCellId(cellId, ijkMin))
604 {
605 bounds[0] = bounds[1] = bounds[2] = bounds[3] = bounds[4] = bounds[5] = 0.0;
606 return;
607 }
608
609 int ijkMax[3];
610 if (!this->GetIJKMaxForIJKMin(ijkMin, ijkMax))
611 {
612 bounds[0] = bounds[1] = bounds[2] = bounds[3] = bounds[4] = bounds[5] = 0.0;
613 return;
614 }
615
616 int loc[3], i, j, k;
617 double xyz[3];
618 const int* extent = this->Extent;
619
620 // Compute the bounds
621 if (ijkMax[2] >= ijkMin[2] && ijkMax[1] >= ijkMin[1] && ijkMax[0] >= ijkMin[0])
622 {
623 bounds[0] = bounds[2] = bounds[4] = VTK_DOUBLE_MAX;
624 bounds[1] = bounds[3] = bounds[5] = VTK_DOUBLE_MIN;
625
626 for (loc[2] = ijkMin[2]; loc[2] <= ijkMax[2]; loc[2]++)
627 {
628 k = loc[2] + extent[4];
629 for (loc[1] = ijkMin[1]; loc[1] <= ijkMax[1]; loc[1]++)
630 {
631 j = loc[1] + extent[2];
632 for (loc[0] = ijkMin[0]; loc[0] <= ijkMax[0]; loc[0]++)
633 {
634 i = loc[0] + extent[0];
635 this->TransformIndexToPhysicalPoint(i, j, k, xyz);
636
637 bounds[0] = (xyz[0] < bounds[0] ? xyz[0] : bounds[0]);
638 bounds[1] = (xyz[0] > bounds[1] ? xyz[0] : bounds[1]);
639 bounds[2] = (xyz[1] < bounds[2] ? xyz[1] : bounds[2]);
640 bounds[3] = (xyz[1] > bounds[3] ? xyz[1] : bounds[3]);
641 bounds[4] = (xyz[2] < bounds[4] ? xyz[2] : bounds[4]);
642 bounds[5] = (xyz[2] > bounds[5] ? xyz[2] : bounds[5]);
643 }
644 }
645 }
646 }
647 else
648 {
649 vtkMath::UninitializeBounds(bounds);
650 }
651 }
652
653 //------------------------------------------------------------------------------
GetPoint(vtkIdType ptId,double x[3])654 void vtkImageData::GetPoint(vtkIdType ptId, double x[3])
655 {
656 const int* extent = this->Extent;
657
658 vtkIdType dims[3];
659 this->GetDimensions(dims);
660
661 x[0] = x[1] = x[2] = 0.0;
662 if (dims[0] == 0 || dims[1] == 0 || dims[2] == 0)
663 {
664 vtkErrorMacro("Requesting a point from an empty image.");
665 return;
666 }
667
668 // "loc" holds the point x,y,z indices
669 int loc[3];
670 loc[0] = loc[1] = loc[2] = 0;
671
672 switch (this->DataDescription)
673 {
674 case VTK_EMPTY:
675 return;
676
677 case VTK_SINGLE_POINT:
678 break;
679
680 case VTK_X_LINE:
681 loc[0] = ptId;
682 break;
683
684 case VTK_Y_LINE:
685 loc[1] = ptId;
686 break;
687
688 case VTK_Z_LINE:
689 loc[2] = ptId;
690 break;
691
692 case VTK_XY_PLANE:
693 loc[0] = ptId % dims[0];
694 loc[1] = ptId / dims[0];
695 break;
696
697 case VTK_YZ_PLANE:
698 loc[1] = ptId % dims[1];
699 loc[2] = ptId / dims[1];
700 break;
701
702 case VTK_XZ_PLANE:
703 loc[0] = ptId % dims[0];
704 loc[2] = ptId / dims[0];
705 break;
706
707 case VTK_XYZ_GRID:
708 loc[0] = ptId % dims[0];
709 loc[1] = (ptId / dims[0]) % dims[1];
710 loc[2] = ptId / (dims[0] * dims[1]);
711 break;
712 }
713
714 int i, j, k;
715 i = loc[0] + extent[0];
716 j = loc[1] + extent[2];
717 k = loc[2] + extent[4];
718 this->TransformIndexToPhysicalPoint(i, j, k, x);
719 }
720
721 //------------------------------------------------------------------------------
FindPoint(double x[3])722 vtkIdType vtkImageData::FindPoint(double x[3])
723 {
724 //
725 // Ensure valid spacing
726 //
727 const double* spacing = this->Spacing;
728 vtkIdType dims[3];
729 this->GetDimensions(dims);
730 std::string ijkLabels[3] = { "I", "J", "K" };
731 for (int i = 0; i < 3; i++)
732 {
733 if (spacing[i] == 0.0 && dims[i] > 1)
734 {
735 vtkWarningMacro("Spacing along the " << ijkLabels[i] << " axis is 0.");
736 return -1;
737 }
738 }
739
740 //
741 // Compute the ijk location
742 //
743 const int* extent = this->Extent;
744 int loc[3];
745 double ijk[3];
746 this->TransformPhysicalPointToContinuousIndex(x, ijk);
747 loc[0] = vtkMath::Floor(ijk[0] + 0.5);
748 loc[1] = vtkMath::Floor(ijk[1] + 0.5);
749 loc[2] = vtkMath::Floor(ijk[2] + 0.5);
750 if (loc[0] < extent[0] || loc[0] > extent[1] || loc[1] < extent[2] || loc[1] > extent[3] ||
751 loc[2] < extent[4] || loc[2] > extent[5])
752 {
753 return -1;
754 }
755 // since point id is relative to the first point actually stored
756 loc[0] -= extent[0];
757 loc[1] -= extent[2];
758 loc[2] -= extent[4];
759
760 //
761 // From this location get the point id
762 //
763 return loc[2] * dims[0] * dims[1] + loc[1] * dims[0] + loc[0];
764 }
765
766 //------------------------------------------------------------------------------
FindCell(double x[3],vtkCell * vtkNotUsed (cell),vtkGenericCell * vtkNotUsed (gencell),vtkIdType vtkNotUsed (cellId),double tol2,int & subId,double pcoords[3],double * weights)767 vtkIdType vtkImageData::FindCell(double x[3], vtkCell* vtkNotUsed(cell),
768 vtkGenericCell* vtkNotUsed(gencell), vtkIdType vtkNotUsed(cellId), double tol2, int& subId,
769 double pcoords[3], double* weights)
770 {
771 return this->FindCell(x, nullptr, 0, tol2, subId, pcoords, weights);
772 }
773
774 //------------------------------------------------------------------------------
FindCell(double x[3],vtkCell * vtkNotUsed (cell),vtkIdType vtkNotUsed (cellId),double tol2,int & subId,double pcoords[3],double * weights)775 vtkIdType vtkImageData::FindCell(double x[3], vtkCell* vtkNotUsed(cell),
776 vtkIdType vtkNotUsed(cellId), double tol2, int& subId, double pcoords[3], double* weights)
777 {
778 int idx[3];
779
780 // Compute the voxel index
781 if (this->ComputeStructuredCoordinates(x, idx, pcoords) == 0)
782 {
783 // If voxel index is out of bounds, check point "x" against the
784 // bounds to see if within tolerance of the bounds.
785 const int* extent = this->Extent;
786 const double* spacing = this->Spacing;
787
788 // Compute squared distance of point x from the boundary
789 double dist2 = 0.0;
790
791 for (int i = 0; i < 3; i++)
792 {
793 int minIdx = extent[i * 2];
794 int maxIdx = extent[i * 2 + 1];
795
796 if (idx[i] < minIdx)
797 {
798 double dist = (idx[i] + pcoords[i] - minIdx) * spacing[i];
799 idx[i] = minIdx;
800 pcoords[i] = 0.0;
801 dist2 += dist * dist;
802 }
803 else if (idx[i] >= maxIdx)
804 {
805 double dist = (idx[i] + pcoords[i] - maxIdx) * spacing[i];
806 if (maxIdx == minIdx)
807 {
808 idx[i] = minIdx;
809 pcoords[i] = 0.0;
810 }
811 else
812 {
813 idx[i] = maxIdx - 1;
814 pcoords[i] = 1.0;
815 }
816 dist2 += dist * dist;
817 }
818 }
819
820 // Check squared distance against the tolerance
821 if (dist2 > tol2)
822 {
823 return -1;
824 }
825 }
826
827 if (weights)
828 {
829 // Shift parametric coordinates for XZ/YZ planes
830 if (this->DataDescription == VTK_XZ_PLANE)
831 {
832 pcoords[1] = pcoords[2];
833 pcoords[2] = 0.0;
834 }
835 else if (this->DataDescription == VTK_YZ_PLANE)
836 {
837 pcoords[0] = pcoords[1];
838 pcoords[1] = pcoords[2];
839 pcoords[2] = 0.0;
840 }
841 else if (this->DataDescription == VTK_XY_PLANE)
842 {
843 pcoords[2] = 0.0;
844 }
845 vtkVoxel::InterpolationFunctions(pcoords, weights);
846 }
847
848 //
849 // From this location get the cell id
850 //
851 subId = 0;
852 return this->ComputeCellId(idx);
853 }
854
855 //------------------------------------------------------------------------------
FindAndGetCell(double x[3],vtkCell * vtkNotUsed (cell),vtkIdType vtkNotUsed (cellId),double tol2,int & subId,double pcoords[3],double * weights)856 vtkCell* vtkImageData::FindAndGetCell(double x[3], vtkCell* vtkNotUsed(cell),
857 vtkIdType vtkNotUsed(cellId), double tol2, int& subId, double pcoords[3], double* weights)
858 {
859 vtkIdType cellId = this->FindCell(x, nullptr, 0, tol2, subId, pcoords, nullptr);
860
861 if (cellId < 0)
862 {
863 return nullptr;
864 }
865
866 vtkCell* cell = this->GetCell(cellId);
867 cell->InterpolateFunctions(pcoords, weights);
868
869 return cell;
870 }
871
872 //------------------------------------------------------------------------------
GetCellType(vtkIdType vtkNotUsed (cellId))873 int vtkImageData::GetCellType(vtkIdType vtkNotUsed(cellId))
874 {
875 switch (this->DataDescription)
876 {
877 case VTK_EMPTY:
878 return VTK_EMPTY_CELL;
879
880 case VTK_SINGLE_POINT:
881 return VTK_VERTEX;
882
883 case VTK_X_LINE:
884 case VTK_Y_LINE:
885 case VTK_Z_LINE:
886 return VTK_LINE;
887
888 case VTK_XY_PLANE:
889 case VTK_YZ_PLANE:
890 case VTK_XZ_PLANE:
891 return VTK_PIXEL;
892
893 case VTK_XYZ_GRID:
894 return VTK_VOXEL;
895
896 default:
897 vtkErrorMacro(<< "Bad data description!");
898 return VTK_EMPTY_CELL;
899 }
900 }
901
902 //------------------------------------------------------------------------------
ComputeBounds()903 void vtkImageData::ComputeBounds()
904 {
905 if (this->GetMTime() <= this->ComputeTime)
906 {
907 return;
908 }
909 const int* extent = this->Extent;
910
911 if (extent[0] > extent[1] || extent[2] > extent[3] || extent[4] > extent[5])
912 {
913 vtkMath::UninitializeBounds(this->Bounds);
914 }
915 else
916 {
917 if (this->DirectionMatrix->IsIdentity())
918 {
919 // Direction is identity: bounds are easy to compute
920 // with only origin and spacing
921 const double* origin = this->Origin;
922 const double* spacing = this->Spacing;
923 int swapXBounds = (spacing[0] < 0); // 1 if true, 0 if false
924 int swapYBounds = (spacing[1] < 0); // 1 if true, 0 if false
925 int swapZBounds = (spacing[2] < 0); // 1 if true, 0 if false
926
927 this->Bounds[0] = origin[0] + (extent[0 + swapXBounds] * spacing[0]);
928 this->Bounds[2] = origin[1] + (extent[2 + swapYBounds] * spacing[1]);
929 this->Bounds[4] = origin[2] + (extent[4 + swapZBounds] * spacing[2]);
930
931 this->Bounds[1] = origin[0] + (extent[1 - swapXBounds] * spacing[0]);
932 this->Bounds[3] = origin[1] + (extent[3 - swapYBounds] * spacing[1]);
933 this->Bounds[5] = origin[2] + (extent[5 - swapZBounds] * spacing[2]);
934 }
935 else
936 {
937 // Direction isn't identity: use IndexToPhysical matrix
938 // to determine the position of the dataset corners
939 int iMin, iMax, jMin, jMax, kMin, kMax;
940 iMin = extent[0];
941 iMax = extent[1];
942 jMin = extent[2];
943 jMax = extent[3];
944 kMin = extent[4];
945 kMax = extent[5];
946 int ijkCorners[8][3] = {
947 { iMin, jMin, kMin },
948 { iMax, jMin, kMin },
949 { iMin, jMax, kMin },
950 { iMax, jMax, kMin },
951 { iMin, jMin, kMax },
952 { iMax, jMin, kMax },
953 { iMin, jMax, kMax },
954 { iMax, jMax, kMax },
955 };
956
957 double xyz[3];
958 double xMin, xMax, yMin, yMax, zMin, zMax;
959 xMin = yMin = zMin = VTK_DOUBLE_MAX;
960 xMax = yMax = zMax = VTK_DOUBLE_MIN;
961 for (int* ijkCorner : ijkCorners)
962 {
963 this->TransformIndexToPhysicalPoint(ijkCorner, xyz);
964 if (xyz[0] < xMin)
965 xMin = xyz[0];
966 if (xyz[0] > xMax)
967 xMax = xyz[0];
968 if (xyz[1] < yMin)
969 yMin = xyz[1];
970 if (xyz[1] > yMax)
971 yMax = xyz[1];
972 if (xyz[2] < zMin)
973 zMin = xyz[2];
974 if (xyz[2] > zMax)
975 zMax = xyz[2];
976 }
977 this->Bounds[0] = xMin;
978 this->Bounds[1] = xMax;
979 this->Bounds[2] = yMin;
980 this->Bounds[3] = yMax;
981 this->Bounds[4] = zMin;
982 this->Bounds[5] = zMax;
983 }
984 }
985 this->ComputeTime.Modified();
986 }
987
988 //------------------------------------------------------------------------------
989 // Given structured coordinates (i,j,k) for a voxel cell, compute the eight
990 // gradient values for the voxel corners. The order in which the gradient
991 // vectors are arranged corresponds to the ordering of the voxel points.
992 // Gradient vector is computed by central differences (except on edges of
993 // volume where forward difference is used). The scalars s are the scalars
994 // from which the gradient is to be computed. This method will treat
995 // only 3D structured point datasets (i.e., volumes).
GetVoxelGradient(int i,int j,int k,vtkDataArray * s,vtkDataArray * g)996 void vtkImageData::GetVoxelGradient(int i, int j, int k, vtkDataArray* s, vtkDataArray* g)
997 {
998 double gv[3];
999 int ii, jj, kk, idx = 0;
1000
1001 for (kk = 0; kk < 2; kk++)
1002 {
1003 for (jj = 0; jj < 2; jj++)
1004 {
1005 for (ii = 0; ii < 2; ii++)
1006 {
1007 this->GetPointGradient(i + ii, j + jj, k + kk, s, gv);
1008 g->SetTuple(idx++, gv);
1009 }
1010 }
1011 }
1012 }
1013
1014 //------------------------------------------------------------------------------
1015 // Given structured coordinates (i,j,k) for a point in a structured point
1016 // dataset, compute the gradient vector from the scalar data at that point.
1017 // The scalars s are the scalars from which the gradient is to be computed.
1018 // This method will treat structured point datasets of any dimension.
GetPointGradient(int i,int j,int k,vtkDataArray * s,double g[3])1019 void vtkImageData::GetPointGradient(int i, int j, int k, vtkDataArray* s, double g[3])
1020 {
1021 const double* ar = this->Spacing;
1022 double sp, sm;
1023 const int* extent = this->Extent;
1024
1025 vtkIdType dims[3];
1026 this->GetDimensions(dims);
1027 vtkIdType ijsize = dims[0] * dims[1];
1028
1029 // Adjust i,j,k to the start of the extent
1030 i -= extent[0];
1031 j -= extent[2];
1032 k -= extent[4];
1033
1034 // Check for out-of-bounds
1035 if (i < 0 || i >= dims[0] || j < 0 || j >= dims[1] || k < 0 || k >= dims[2])
1036 {
1037 g[0] = g[1] = g[2] = 0.0;
1038 return;
1039 }
1040
1041 // i-axis
1042 if (dims[0] == 1)
1043 {
1044 g[0] = 0.0;
1045 }
1046 else if (i == 0)
1047 {
1048 sp = s->GetComponent(i + 1 + j * dims[0] + k * ijsize, 0);
1049 sm = s->GetComponent(i + j * dims[0] + k * ijsize, 0);
1050 g[0] = (sm - sp) / ar[0];
1051 }
1052 else if (i == (dims[0] - 1))
1053 {
1054 sp = s->GetComponent(i + j * dims[0] + k * ijsize, 0);
1055 sm = s->GetComponent(i - 1 + j * dims[0] + k * ijsize, 0);
1056 g[0] = (sm - sp) / ar[0];
1057 }
1058 else
1059 {
1060 sp = s->GetComponent(i + 1 + j * dims[0] + k * ijsize, 0);
1061 sm = s->GetComponent(i - 1 + j * dims[0] + k * ijsize, 0);
1062 g[0] = 0.5 * (sm - sp) / ar[0];
1063 }
1064
1065 // j-axis
1066 if (dims[1] == 1)
1067 {
1068 g[1] = 0.0;
1069 }
1070 else if (j == 0)
1071 {
1072 sp = s->GetComponent(i + (j + 1) * dims[0] + k * ijsize, 0);
1073 sm = s->GetComponent(i + j * dims[0] + k * ijsize, 0);
1074 g[1] = (sm - sp) / ar[1];
1075 }
1076 else if (j == (dims[1] - 1))
1077 {
1078 sp = s->GetComponent(i + j * dims[0] + k * ijsize, 0);
1079 sm = s->GetComponent(i + (j - 1) * dims[0] + k * ijsize, 0);
1080 g[1] = (sm - sp) / ar[1];
1081 }
1082 else
1083 {
1084 sp = s->GetComponent(i + (j + 1) * dims[0] + k * ijsize, 0);
1085 sm = s->GetComponent(i + (j - 1) * dims[0] + k * ijsize, 0);
1086 g[1] = 0.5 * (sm - sp) / ar[1];
1087 }
1088
1089 // k-axis
1090 if (dims[2] == 1)
1091 {
1092 g[2] = 0.0;
1093 }
1094 else if (k == 0)
1095 {
1096 sp = s->GetComponent(i + j * dims[0] + (k + 1) * ijsize, 0);
1097 sm = s->GetComponent(i + j * dims[0] + k * ijsize, 0);
1098 g[2] = (sm - sp) / ar[2];
1099 }
1100 else if (k == (dims[2] - 1))
1101 {
1102 sp = s->GetComponent(i + j * dims[0] + k * ijsize, 0);
1103 sm = s->GetComponent(i + j * dims[0] + (k - 1) * ijsize, 0);
1104 g[2] = (sm - sp) / ar[2];
1105 }
1106 else
1107 {
1108 sp = s->GetComponent(i + j * dims[0] + (k + 1) * ijsize, 0);
1109 sm = s->GetComponent(i + j * dims[0] + (k - 1) * ijsize, 0);
1110 g[2] = 0.5 * (sm - sp) / ar[2];
1111 }
1112
1113 // Apply direction transform to get in xyz coordinate system
1114 // Note: we already applied the spacing when handling the ijk
1115 // axis above, and do not need to translate by the origin
1116 // since this is a gradient computation
1117 this->DirectionMatrix->MultiplyPoint(g, g);
1118 }
1119
1120 //------------------------------------------------------------------------------
1121 // Set dimensions of structured points dataset.
SetDimensions(int i,int j,int k)1122 void vtkImageData::SetDimensions(int i, int j, int k)
1123 {
1124 this->SetExtent(0, i - 1, 0, j - 1, 0, k - 1);
1125 }
1126
1127 //------------------------------------------------------------------------------
1128 // Set dimensions of structured points dataset.
SetDimensions(const int dim[3])1129 void vtkImageData::SetDimensions(const int dim[3])
1130 {
1131 this->SetExtent(0, dim[0] - 1, 0, dim[1] - 1, 0, dim[2] - 1);
1132 }
1133
1134 //------------------------------------------------------------------------------
1135 // Convenience function computes the structured coordinates for a point x[3].
1136 // The voxel is specified by the array ijk[3], and the parametric coordinates
1137 // in the cell are specified with pcoords[3]. The function returns a 0 if the
1138 // point x is outside of the volume, and a 1 if inside the volume.
ComputeStructuredCoordinates(const double x[3],int ijk[3],double pcoords[3])1139 int vtkImageData::ComputeStructuredCoordinates(const double x[3], int ijk[3], double pcoords[3])
1140 {
1141 // tolerance is needed for floating points error margin
1142 // (this is squared tolerance)
1143 const double tol2 = 1e-12;
1144
1145 //
1146 // Compute the ijk location
1147 //
1148 double doubleLoc[3];
1149 this->TransformPhysicalPointToContinuousIndex(x, doubleLoc);
1150
1151 const int* extent = this->Extent;
1152 int isInBounds = 1;
1153 for (int i = 0; i < 3; i++)
1154 {
1155 // Floor for negative indexes.
1156 ijk[i] = vtkMath::Floor(doubleLoc[i]); // integer
1157 pcoords[i] = doubleLoc[i] - ijk[i]; // >= 0 and < 1
1158
1159 int tmpInBounds = 0;
1160 int minExt = extent[i * 2];
1161 int maxExt = extent[i * 2 + 1];
1162
1163 // check if data is one pixel thick as well as
1164 // low boundary check
1165 if (minExt == maxExt || ijk[i] < minExt)
1166 {
1167 double dist = doubleLoc[i] - minExt;
1168 if (dist * dist <= tol2)
1169 {
1170 pcoords[i] = 0.0;
1171 ijk[i] = minExt;
1172 tmpInBounds = 1;
1173 }
1174 }
1175
1176 // high boundary check
1177 else if (ijk[i] >= maxExt)
1178 {
1179 double dist = doubleLoc[i] - maxExt;
1180 if (dist * dist <= tol2)
1181 {
1182 // make sure index is within the allowed cell index range
1183 pcoords[i] = 1.0;
1184 ijk[i] = maxExt - 1;
1185 tmpInBounds = 1;
1186 }
1187 }
1188
1189 // else index is definitely within bounds
1190 else
1191 {
1192 tmpInBounds = 1;
1193 }
1194
1195 // clear isInBounds if out of bounds for this dimension
1196 isInBounds = (isInBounds & tmpInBounds);
1197 }
1198
1199 return isInBounds;
1200 }
1201
1202 //------------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)1203 void vtkImageData::PrintSelf(ostream& os, vtkIndent indent)
1204 {
1205 this->Superclass::PrintSelf(os, indent);
1206
1207 int idx;
1208 const double* direction = this->GetDirectionMatrix()->GetData();
1209 const int* dims = this->GetDimensions();
1210 const int* extent = this->Extent;
1211
1212 os << indent << "Spacing: (" << this->Spacing[0] << ", " << this->Spacing[1] << ", "
1213 << this->Spacing[2] << ")\n";
1214 os << indent << "Origin: (" << this->Origin[0] << ", " << this->Origin[1] << ", "
1215 << this->Origin[2] << ")\n";
1216 os << indent << "Direction: (" << direction[0];
1217 for (idx = 1; idx < 9; ++idx)
1218 {
1219 os << ", " << direction[idx];
1220 }
1221 os << ")\n";
1222 os << indent << "Dimensions: (" << dims[0] << ", " << dims[1] << ", " << dims[2] << ")\n";
1223 os << indent << "Increments: (" << this->Increments[0] << ", " << this->Increments[1] << ", "
1224 << this->Increments[2] << ")\n";
1225 os << indent << "Extent: (" << extent[0];
1226 for (idx = 1; idx < 6; ++idx)
1227 {
1228 os << ", " << extent[idx];
1229 }
1230 os << ")\n";
1231 }
1232
1233 //------------------------------------------------------------------------------
SetNumberOfScalarComponents(int num,vtkInformation * meta_data)1234 void vtkImageData::SetNumberOfScalarComponents(int num, vtkInformation* meta_data)
1235 {
1236 vtkDataObject::SetPointDataActiveScalarInfo(meta_data, -1, num);
1237 }
1238
1239 //------------------------------------------------------------------------------
HasNumberOfScalarComponents(vtkInformation * meta_data)1240 bool vtkImageData::HasNumberOfScalarComponents(vtkInformation* meta_data)
1241 {
1242 vtkInformation* scalarInfo = vtkDataObject::GetActiveFieldInformation(
1243 meta_data, FIELD_ASSOCIATION_POINTS, vtkDataSetAttributes::SCALARS);
1244 if (!scalarInfo)
1245 {
1246 return false;
1247 }
1248 return scalarInfo->Has(FIELD_NUMBER_OF_COMPONENTS()) != 0;
1249 }
1250
1251 //------------------------------------------------------------------------------
GetNumberOfScalarComponents(vtkInformation * meta_data)1252 int vtkImageData::GetNumberOfScalarComponents(vtkInformation* meta_data)
1253 {
1254 vtkInformation* scalarInfo = vtkDataObject::GetActiveFieldInformation(
1255 meta_data, FIELD_ASSOCIATION_POINTS, vtkDataSetAttributes::SCALARS);
1256 if (scalarInfo && scalarInfo->Has(FIELD_NUMBER_OF_COMPONENTS()))
1257 {
1258 return scalarInfo->Get(FIELD_NUMBER_OF_COMPONENTS());
1259 }
1260 return 1;
1261 }
1262
1263 //------------------------------------------------------------------------------
GetNumberOfScalarComponents()1264 int vtkImageData::GetNumberOfScalarComponents()
1265 {
1266 vtkDataArray* scalars = this->GetPointData()->GetScalars();
1267 if (scalars)
1268 {
1269 return scalars->GetNumberOfComponents();
1270 }
1271 return 1;
1272 }
1273
1274 //------------------------------------------------------------------------------
GetIncrements()1275 vtkIdType* vtkImageData::GetIncrements()
1276 {
1277 // Make sure the increments are up to date. The filter bypass and update
1278 // mechanism make it tricky to update the increments anywhere other than here
1279 this->ComputeIncrements();
1280
1281 return this->Increments;
1282 }
1283
1284 //------------------------------------------------------------------------------
GetIncrements(vtkDataArray * scalars)1285 vtkIdType* vtkImageData::GetIncrements(vtkDataArray* scalars)
1286 {
1287 // Make sure the increments are up to date. The filter bypass and update
1288 // mechanism make it tricky to update the increments anywhere other than here
1289 this->ComputeIncrements(scalars);
1290
1291 return this->Increments;
1292 }
1293
1294 //------------------------------------------------------------------------------
GetIncrements(vtkIdType & incX,vtkIdType & incY,vtkIdType & incZ)1295 void vtkImageData::GetIncrements(vtkIdType& incX, vtkIdType& incY, vtkIdType& incZ)
1296 {
1297 vtkIdType inc[3];
1298 this->ComputeIncrements(inc);
1299 incX = inc[0];
1300 incY = inc[1];
1301 incZ = inc[2];
1302 }
1303
1304 //------------------------------------------------------------------------------
GetIncrements(vtkDataArray * scalars,vtkIdType & incX,vtkIdType & incY,vtkIdType & incZ)1305 void vtkImageData::GetIncrements(
1306 vtkDataArray* scalars, vtkIdType& incX, vtkIdType& incY, vtkIdType& incZ)
1307 {
1308 vtkIdType inc[3];
1309 this->ComputeIncrements(scalars, inc);
1310 incX = inc[0];
1311 incY = inc[1];
1312 incZ = inc[2];
1313 }
1314
1315 //------------------------------------------------------------------------------
GetIncrements(vtkIdType inc[3])1316 void vtkImageData::GetIncrements(vtkIdType inc[3])
1317 {
1318 this->ComputeIncrements(inc);
1319 }
1320
1321 //------------------------------------------------------------------------------
GetIncrements(vtkDataArray * scalars,vtkIdType inc[3])1322 void vtkImageData::GetIncrements(vtkDataArray* scalars, vtkIdType inc[3])
1323 {
1324 this->ComputeIncrements(scalars, inc);
1325 }
1326
1327 //------------------------------------------------------------------------------
GetContinuousIncrements(int extent[6],vtkIdType & incX,vtkIdType & incY,vtkIdType & incZ)1328 void vtkImageData::GetContinuousIncrements(
1329 int extent[6], vtkIdType& incX, vtkIdType& incY, vtkIdType& incZ)
1330 {
1331 this->GetContinuousIncrements(this->GetPointData()->GetScalars(), extent, incX, incY, incZ);
1332 }
1333 //------------------------------------------------------------------------------
GetContinuousIncrements(vtkDataArray * scalars,int extent[6],vtkIdType & incX,vtkIdType & incY,vtkIdType & incZ)1334 void vtkImageData::GetContinuousIncrements(
1335 vtkDataArray* scalars, int extent[6], vtkIdType& incX, vtkIdType& incY, vtkIdType& incZ)
1336 {
1337 int e0, e1, e2, e3;
1338
1339 incX = 0;
1340 const int* selfExtent = this->Extent;
1341
1342 e0 = extent[0];
1343 if (e0 < selfExtent[0])
1344 {
1345 e0 = selfExtent[0];
1346 }
1347 e1 = extent[1];
1348 if (e1 > selfExtent[1])
1349 {
1350 e1 = selfExtent[1];
1351 }
1352 e2 = extent[2];
1353 if (e2 < selfExtent[2])
1354 {
1355 e2 = selfExtent[2];
1356 }
1357 e3 = extent[3];
1358 if (e3 > selfExtent[3])
1359 {
1360 e3 = selfExtent[3];
1361 }
1362
1363 // Make sure the increments are up to date
1364 vtkIdType inc[3];
1365 this->ComputeIncrements(scalars, inc);
1366
1367 incY = inc[1] - (e1 - e0 + 1) * inc[0];
1368 incZ = inc[2] - (e3 - e2 + 1) * inc[1];
1369 }
1370
1371 //------------------------------------------------------------------------------
1372 // This method computes the increments from the MemoryOrder and the extent.
1373 // This version assumes we are using the Active Scalars
ComputeIncrements(vtkIdType inc[3])1374 void vtkImageData::ComputeIncrements(vtkIdType inc[3])
1375 {
1376 this->ComputeIncrements(this->GetPointData()->GetScalars(), inc);
1377 }
1378
1379 //------------------------------------------------------------------------------
1380 // This method computes the increments from the MemoryOrder and the extent.
ComputeIncrements(vtkDataArray * scalars,vtkIdType inc[3])1381 void vtkImageData::ComputeIncrements(vtkDataArray* scalars, vtkIdType inc[3])
1382 {
1383 if (!scalars)
1384 {
1385 vtkErrorMacro("No Scalar Field has been specified - assuming 1 component!");
1386 this->ComputeIncrements(1, inc);
1387 }
1388 else
1389 {
1390 this->ComputeIncrements(scalars->GetNumberOfComponents(), inc);
1391 }
1392 }
1393 //------------------------------------------------------------------------------
1394 // This method computes the increments from the MemoryOrder and the extent.
ComputeIncrements(int numberOfComponents,vtkIdType inc[3])1395 void vtkImageData::ComputeIncrements(int numberOfComponents, vtkIdType inc[3])
1396 {
1397 int idx;
1398 vtkIdType incr = numberOfComponents;
1399 const int* extent = this->Extent;
1400
1401 for (idx = 0; idx < 3; ++idx)
1402 {
1403 inc[idx] = incr;
1404 incr *= (extent[idx * 2 + 1] - extent[idx * 2] + 1);
1405 }
1406 }
1407
1408 //------------------------------------------------------------------------------
1409 template <class TIn, class TOut>
vtkImageDataConvertScalar(TIn * in,TOut * out)1410 void vtkImageDataConvertScalar(TIn* in, TOut* out)
1411 {
1412 *out = static_cast<TOut>(*in);
1413 }
1414
1415 //------------------------------------------------------------------------------
GetScalarComponentAsDouble(int x,int y,int z,int comp)1416 double vtkImageData::GetScalarComponentAsDouble(int x, int y, int z, int comp)
1417 {
1418 // Check the component index.
1419 if (comp < 0 || comp >= this->GetNumberOfScalarComponents())
1420 {
1421 vtkErrorMacro("Bad component index " << comp);
1422 return 0.0;
1423 }
1424
1425 vtkIdType index = this->GetScalarIndex(x, y, z);
1426 if (index < 0)
1427 {
1428 // An error message was already generated by GetScalarIndex.
1429 return 0.0;
1430 }
1431
1432 vtkDataArray* scalars = this->GetPointData()->GetScalars();
1433 return scalars->GetComponent(index, comp);
1434 }
1435
1436 //------------------------------------------------------------------------------
SetScalarComponentFromDouble(int x,int y,int z,int comp,double value)1437 void vtkImageData::SetScalarComponentFromDouble(int x, int y, int z, int comp, double value)
1438 {
1439 // Check the component index.
1440 if (comp < 0 || comp >= this->GetNumberOfScalarComponents())
1441 {
1442 vtkErrorMacro("Bad component index " << comp);
1443 return;
1444 }
1445
1446 vtkIdType index = this->GetScalarIndex(x, y, z);
1447 if (index < 0)
1448 {
1449 // An error message was already generated by GetScalarIndex.
1450 return;
1451 }
1452
1453 vtkDataArray* scalars = this->GetPointData()->GetScalars();
1454 scalars->SetComponent(index, comp, value);
1455 }
1456
1457 //------------------------------------------------------------------------------
GetScalarComponentAsFloat(int x,int y,int z,int comp)1458 float vtkImageData::GetScalarComponentAsFloat(int x, int y, int z, int comp)
1459 {
1460 return this->GetScalarComponentAsDouble(x, y, z, comp);
1461 }
1462
1463 //------------------------------------------------------------------------------
SetScalarComponentFromFloat(int x,int y,int z,int comp,float value)1464 void vtkImageData::SetScalarComponentFromFloat(int x, int y, int z, int comp, float value)
1465 {
1466 this->SetScalarComponentFromDouble(x, y, z, comp, value);
1467 }
1468
1469 //------------------------------------------------------------------------------
1470 // This Method returns a pointer to a location in the vtkImageData.
1471 // Coordinates are in pixel units and are relative to the whole
1472 // image origin.
GetScalarPointer(int x,int y,int z)1473 void* vtkImageData::GetScalarPointer(int x, int y, int z)
1474 {
1475 int tmp[3];
1476 tmp[0] = x;
1477 tmp[1] = y;
1478 tmp[2] = z;
1479 return this->GetScalarPointer(tmp);
1480 }
1481
1482 //------------------------------------------------------------------------------
1483 // This Method returns a pointer to a location in the vtkImageData.
1484 // Coordinates are in pixel units and are relative to the whole
1485 // image origin.
GetScalarPointerForExtent(int extent[6])1486 void* vtkImageData::GetScalarPointerForExtent(int extent[6])
1487 {
1488 int tmp[3];
1489 tmp[0] = extent[0];
1490 tmp[1] = extent[2];
1491 tmp[2] = extent[4];
1492 return this->GetScalarPointer(tmp);
1493 }
1494
1495 //------------------------------------------------------------------------------
GetScalarPointer(int coordinate[3])1496 void* vtkImageData::GetScalarPointer(int coordinate[3])
1497 {
1498 vtkDataArray* scalars = this->GetPointData()->GetScalars();
1499
1500 // Make sure the array has been allocated.
1501 if (scalars == nullptr)
1502 {
1503 // vtkDebugMacro("Allocating scalars in ImageData");
1504 // abort();
1505 // this->AllocateScalars();
1506 // scalars = this->PointData->GetScalars();
1507 return nullptr;
1508 }
1509
1510 const int* extent = this->Extent;
1511 // error checking: since most access will be from pointer arithmetic.
1512 // this should not waste much time.
1513 for (int idx = 0; idx < 3; ++idx)
1514 {
1515 if (coordinate[idx] < extent[idx * 2] || coordinate[idx] > extent[idx * 2 + 1])
1516 {
1517 vtkErrorMacro(<< "GetScalarPointer: Pixel (" << coordinate[0] << ", " << coordinate[1] << ", "
1518 << coordinate[2] << ") not in memory.\n Current extent= (" << extent[0] << ", "
1519 << extent[1] << ", " << extent[2] << ", " << extent[3] << ", " << extent[4]
1520 << ", " << extent[5] << ")");
1521 return nullptr;
1522 }
1523 }
1524
1525 return this->GetArrayPointer(scalars, coordinate);
1526 }
1527
1528 //------------------------------------------------------------------------------
1529 // This method returns a pointer to the origin of the vtkImageData.
GetScalarPointer()1530 void* vtkImageData::GetScalarPointer()
1531 {
1532 if (this->PointData->GetScalars() == nullptr)
1533 {
1534 // vtkDebugMacro("Allocating scalars in ImageData");
1535 // abort();
1536 // this->AllocateScalars();
1537 return nullptr;
1538 }
1539 return this->PointData->GetScalars()->GetVoidPointer(0);
1540 }
1541
1542 //------------------------------------------------------------------------------
1543 // This Method returns an index to a location in the vtkImageData.
1544 // Coordinates are in pixel units and are relative to the whole
1545 // image origin.
GetScalarIndex(int x,int y,int z)1546 vtkIdType vtkImageData::GetScalarIndex(int x, int y, int z)
1547 {
1548 int tmp[3];
1549 tmp[0] = x;
1550 tmp[1] = y;
1551 tmp[2] = z;
1552 return this->GetScalarIndex(tmp);
1553 }
1554
1555 //------------------------------------------------------------------------------
1556 // This Method returns an index to a location in the vtkImageData.
1557 // Coordinates are in pixel units and are relative to the whole
1558 // image origin.
GetScalarIndexForExtent(int extent[6])1559 vtkIdType vtkImageData::GetScalarIndexForExtent(int extent[6])
1560 {
1561 int tmp[3];
1562 tmp[0] = extent[0];
1563 tmp[1] = extent[2];
1564 tmp[2] = extent[4];
1565 return this->GetScalarIndex(tmp);
1566 }
1567
1568 //------------------------------------------------------------------------------
GetScalarIndex(int coordinate[3])1569 vtkIdType vtkImageData::GetScalarIndex(int coordinate[3])
1570 {
1571 vtkDataArray* scalars = this->GetPointData()->GetScalars();
1572
1573 // Make sure the array has been allocated.
1574 if (scalars == nullptr)
1575 {
1576 return -1;
1577 }
1578
1579 const int* extent = this->Extent;
1580 // error checking: since most access will be from pointer arithmetic.
1581 // this should not waste much time.
1582 for (int idx = 0; idx < 3; ++idx)
1583 {
1584 if (coordinate[idx] < extent[idx * 2] || coordinate[idx] > extent[idx * 2 + 1])
1585 {
1586 vtkErrorMacro(<< "GetScalarIndex: Pixel (" << coordinate[0] << ", " << coordinate[1] << ", "
1587 << coordinate[2] << ") not in memory.\n Current extent= (" << extent[0] << ", "
1588 << extent[1] << ", " << extent[2] << ", " << extent[3] << ", " << extent[4]
1589 << ", " << extent[5] << ")");
1590 return -1;
1591 }
1592 }
1593
1594 return this->GetTupleIndex(scalars, coordinate);
1595 }
1596
1597 //------------------------------------------------------------------------------
SetScalarType(int type,vtkInformation * meta_data)1598 void vtkImageData::SetScalarType(int type, vtkInformation* meta_data)
1599 {
1600 vtkDataObject::SetPointDataActiveScalarInfo(meta_data, type, -1);
1601 }
1602
1603 //------------------------------------------------------------------------------
GetScalarType()1604 int vtkImageData::GetScalarType()
1605 {
1606 vtkDataArray* scalars = this->GetPointData()->GetScalars();
1607 if (!scalars)
1608 {
1609 return VTK_DOUBLE;
1610 }
1611 return scalars->GetDataType();
1612 }
1613
1614 //------------------------------------------------------------------------------
HasScalarType(vtkInformation * meta_data)1615 bool vtkImageData::HasScalarType(vtkInformation* meta_data)
1616 {
1617 vtkInformation* scalarInfo = vtkDataObject::GetActiveFieldInformation(
1618 meta_data, FIELD_ASSOCIATION_POINTS, vtkDataSetAttributes::SCALARS);
1619 if (!scalarInfo)
1620 {
1621 return false;
1622 }
1623
1624 return scalarInfo->Has(FIELD_ARRAY_TYPE()) != 0;
1625 }
1626
1627 //------------------------------------------------------------------------------
GetScalarType(vtkInformation * meta_data)1628 int vtkImageData::GetScalarType(vtkInformation* meta_data)
1629 {
1630 vtkInformation* scalarInfo = vtkDataObject::GetActiveFieldInformation(
1631 meta_data, FIELD_ASSOCIATION_POINTS, vtkDataSetAttributes::SCALARS);
1632 if (scalarInfo)
1633 {
1634 return scalarInfo->Get(FIELD_ARRAY_TYPE());
1635 }
1636 return VTK_DOUBLE;
1637 }
1638
1639 //------------------------------------------------------------------------------
AllocateScalars(vtkInformation * pipeline_info)1640 void vtkImageData::AllocateScalars(vtkInformation* pipeline_info)
1641 {
1642 auto mkhold = vtkMemkindRAII(this->GetIsInMemkind());
1643 int newType = VTK_DOUBLE;
1644 int newNumComp = 1;
1645
1646 if (pipeline_info)
1647 {
1648 vtkInformation* scalarInfo = vtkDataObject::GetActiveFieldInformation(
1649 pipeline_info, FIELD_ASSOCIATION_POINTS, vtkDataSetAttributes::SCALARS);
1650 if (scalarInfo)
1651 {
1652 newType = scalarInfo->Get(FIELD_ARRAY_TYPE());
1653 if (scalarInfo->Has(FIELD_NUMBER_OF_COMPONENTS()))
1654 {
1655 newNumComp = scalarInfo->Get(FIELD_NUMBER_OF_COMPONENTS());
1656 }
1657 }
1658 }
1659
1660 this->AllocateScalars(newType, newNumComp);
1661 }
1662
1663 //------------------------------------------------------------------------------
AllocateScalars(int dataType,int numComponents)1664 void vtkImageData::AllocateScalars(int dataType, int numComponents)
1665 {
1666 auto mkhold = vtkMemkindRAII(this->GetIsInMemkind());
1667 vtkDataArray* scalars;
1668
1669 // if the scalar type has not been set then we have a problem
1670 if (dataType == VTK_VOID)
1671 {
1672 vtkErrorMacro("Attempt to allocate scalars before scalar type was set!.");
1673 return;
1674 }
1675
1676 const int* extent = this->Extent;
1677 // Use vtkIdType to avoid overflow on large images
1678 vtkIdType dims[3];
1679 dims[0] = extent[1] - extent[0] + 1;
1680 dims[1] = extent[3] - extent[2] + 1;
1681 dims[2] = extent[5] - extent[4] + 1;
1682 vtkIdType imageSize = dims[0] * dims[1] * dims[2];
1683
1684 // if we currently have scalars then just adjust the size
1685 scalars = this->PointData->GetScalars();
1686 if (scalars && scalars->GetDataType() == dataType && scalars->GetReferenceCount() == 1)
1687 {
1688 scalars->SetNumberOfComponents(numComponents);
1689 scalars->SetNumberOfTuples(imageSize);
1690 // Since the execute method will be modifying the scalars
1691 // directly.
1692 scalars->Modified();
1693 return;
1694 }
1695
1696 // allocate the new scalars
1697 scalars = vtkDataArray::CreateDataArray(dataType);
1698 scalars->SetNumberOfComponents(numComponents);
1699 scalars->SetName("ImageScalars");
1700
1701 // allocate enough memory
1702 scalars->SetNumberOfTuples(imageSize);
1703
1704 this->PointData->SetScalars(scalars);
1705 scalars->Delete();
1706 }
1707
1708 //------------------------------------------------------------------------------
GetScalarSize(vtkInformation * meta_data)1709 int vtkImageData::GetScalarSize(vtkInformation* meta_data)
1710 {
1711 return vtkDataArray::GetDataTypeSize(this->GetScalarType(meta_data));
1712 }
1713
GetScalarSize()1714 int vtkImageData::GetScalarSize()
1715 {
1716 vtkDataArray* scalars = this->GetPointData()->GetScalars();
1717 if (!scalars)
1718 {
1719 return vtkDataArray::GetDataTypeSize(VTK_DOUBLE);
1720 }
1721 return vtkDataArray::GetDataTypeSize(scalars->GetDataType());
1722 }
1723
1724 //------------------------------------------------------------------------------
1725 // This templated function executes the filter for any type of data.
1726 template <class IT, class OT>
vtkImageDataCastExecute(vtkImageData * inData,IT * inPtr,vtkImageData * outData,OT * outPtr,int outExt[6])1727 void vtkImageDataCastExecute(
1728 vtkImageData* inData, IT* inPtr, vtkImageData* outData, OT* outPtr, int outExt[6])
1729 {
1730 int idxR, idxY, idxZ;
1731 int maxY, maxZ;
1732 vtkIdType inIncX, inIncY, inIncZ;
1733 vtkIdType outIncX, outIncY, outIncZ;
1734 int rowLength;
1735
1736 // find the region to loop over
1737 rowLength = (outExt[1] - outExt[0] + 1) * inData->GetNumberOfScalarComponents();
1738 maxY = outExt[3] - outExt[2];
1739 maxZ = outExt[5] - outExt[4];
1740
1741 // Get increments to march through data
1742 inData->GetContinuousIncrements(outExt, inIncX, inIncY, inIncZ);
1743 outData->GetContinuousIncrements(outExt, outIncX, outIncY, outIncZ);
1744
1745 // Loop through output pixels
1746 for (idxZ = 0; idxZ <= maxZ; idxZ++)
1747 {
1748 for (idxY = 0; idxY <= maxY; idxY++)
1749 {
1750 for (idxR = 0; idxR < rowLength; idxR++)
1751 {
1752 // Pixel operation
1753 *outPtr = static_cast<OT>(*inPtr);
1754 outPtr++;
1755 inPtr++;
1756 }
1757 outPtr += outIncY;
1758 inPtr += inIncY;
1759 }
1760 outPtr += outIncZ;
1761 inPtr += inIncZ;
1762 }
1763 }
1764
1765 //------------------------------------------------------------------------------
1766 template <class T>
vtkImageDataCastExecute(vtkImageData * inData,T * inPtr,vtkImageData * outData,int outExt[6])1767 void vtkImageDataCastExecute(vtkImageData* inData, T* inPtr, vtkImageData* outData, int outExt[6])
1768 {
1769 void* outPtr = outData->GetScalarPointerForExtent(outExt);
1770
1771 if (outPtr == nullptr)
1772 {
1773 vtkGenericWarningMacro("Scalars not allocated.");
1774 return;
1775 }
1776
1777 int scalarType = outData->GetPointData()->GetScalars()->GetDataType();
1778 switch (scalarType)
1779 {
1780 vtkTemplateMacro(vtkImageDataCastExecute(
1781 inData, static_cast<T*>(inPtr), outData, static_cast<VTK_TT*>(outPtr), outExt));
1782 default:
1783 vtkGenericWarningMacro("Execute: Unknown output ScalarType");
1784 return;
1785 }
1786 }
1787
1788 //------------------------------------------------------------------------------
1789 // This method is passed a input and output region, and executes the filter
1790 // algorithm to fill the output from the input.
1791 // It just executes a switch statement to call the correct function for
1792 // the regions data types.
CopyAndCastFrom(vtkImageData * inData,int extent[6])1793 void vtkImageData::CopyAndCastFrom(vtkImageData* inData, int extent[6])
1794 {
1795 void* inPtr = inData->GetScalarPointerForExtent(extent);
1796
1797 if (inPtr == nullptr)
1798 {
1799 vtkErrorMacro("Scalars not allocated.");
1800 return;
1801 }
1802
1803 int scalarType = inData->GetPointData()->GetScalars()->GetDataType();
1804 switch (scalarType)
1805 {
1806 vtkTemplateMacro(vtkImageDataCastExecute(inData, static_cast<VTK_TT*>(inPtr), this, extent));
1807 default:
1808 vtkErrorMacro(<< "Execute: Unknown input ScalarType");
1809 return;
1810 }
1811 }
1812
1813 //------------------------------------------------------------------------------
Crop(const int * updateExtent)1814 void vtkImageData::Crop(const int* updateExtent)
1815 {
1816 // Do nothing for empty datasets:
1817 for (int dim = 0; dim < 3; ++dim)
1818 {
1819 if (this->Extent[2 * dim] > this->Extent[2 * dim + 1])
1820 {
1821 vtkDebugMacro(<< "Refusing to crop empty dataset.");
1822 return;
1823 }
1824 }
1825
1826 int nExt[6];
1827 int idxX, idxY, idxZ;
1828 int maxX, maxY, maxZ;
1829 vtkIdType outId, inId, inIdY, inIdZ, incZ, incY;
1830 vtkImageData* newImage;
1831 vtkIdType numPts, numCells, tmp;
1832 const int* extent = this->Extent;
1833
1834 // If extents already match, then we need to do nothing.
1835 if (extent[0] == updateExtent[0] && extent[1] == updateExtent[1] &&
1836 extent[2] == updateExtent[2] && extent[3] == updateExtent[3] && extent[4] == updateExtent[4] &&
1837 extent[5] == updateExtent[5])
1838 {
1839 return;
1840 }
1841
1842 // Take the intersection of the two extent so that
1843 // we are not asking for more than the extent.
1844 memcpy(nExt, updateExtent, 6 * sizeof(int));
1845 if (nExt[0] < extent[0])
1846 {
1847 nExt[0] = extent[0];
1848 }
1849 if (nExt[1] > extent[1])
1850 {
1851 nExt[1] = extent[1];
1852 }
1853 if (nExt[2] < extent[2])
1854 {
1855 nExt[2] = extent[2];
1856 }
1857 if (nExt[3] > extent[3])
1858 {
1859 nExt[3] = extent[3];
1860 }
1861 if (nExt[4] < extent[4])
1862 {
1863 nExt[4] = extent[4];
1864 }
1865 if (nExt[5] > extent[5])
1866 {
1867 nExt[5] = extent[5];
1868 }
1869
1870 // If the extents are the same just return.
1871 if (extent[0] == nExt[0] && extent[1] == nExt[1] && extent[2] == nExt[2] &&
1872 extent[3] == nExt[3] && extent[4] == nExt[4] && extent[5] == nExt[5])
1873 {
1874 vtkDebugMacro("Extents already match.");
1875 return;
1876 }
1877
1878 // How many point/cells.
1879 numPts = (nExt[1] - nExt[0] + 1) * (nExt[3] - nExt[2] + 1) * (nExt[5] - nExt[4] + 1);
1880 // Conditional are to handle 3d, 2d, and even 1d images.
1881 tmp = nExt[1] - nExt[0];
1882 if (tmp <= 0)
1883 {
1884 tmp = 1;
1885 }
1886 numCells = tmp;
1887 tmp = nExt[3] - nExt[2];
1888 if (tmp <= 0)
1889 {
1890 tmp = 1;
1891 }
1892 numCells *= tmp;
1893 tmp = nExt[5] - nExt[4];
1894 if (tmp <= 0)
1895 {
1896 tmp = 1;
1897 }
1898 numCells *= tmp;
1899
1900 // Create a new temporary image.
1901 newImage = vtkImageData::New();
1902 newImage->SetExtent(nExt);
1903 vtkPointData* npd = newImage->GetPointData();
1904 vtkCellData* ncd = newImage->GetCellData();
1905 npd->CopyAllocate(this->PointData, numPts);
1906 ncd->CopyAllocate(this->CellData, numCells);
1907
1908 // Loop through outData points
1909 incY = extent[1] - extent[0] + 1;
1910 incZ = (extent[3] - extent[2] + 1) * incY;
1911 outId = 0;
1912 inIdZ = incZ * (nExt[4] - extent[4]) + incY * (nExt[2] - extent[2]) + (nExt[0] - extent[0]);
1913
1914 for (idxZ = nExt[4]; idxZ <= nExt[5]; idxZ++)
1915 {
1916 inIdY = inIdZ;
1917 for (idxY = nExt[2]; idxY <= nExt[3]; idxY++)
1918 {
1919 inId = inIdY;
1920 for (idxX = nExt[0]; idxX <= nExt[1]; idxX++)
1921 {
1922 npd->CopyData(this->PointData, inId, outId);
1923 ++inId;
1924 ++outId;
1925 }
1926 inIdY += incY;
1927 }
1928 inIdZ += incZ;
1929 }
1930
1931 // Loop through outData cells
1932 // Have to handle the 2d and 1d cases.
1933 maxX = nExt[1];
1934 maxY = nExt[3];
1935 maxZ = nExt[5];
1936 if (maxX == nExt[0])
1937 {
1938 ++maxX;
1939 }
1940 if (maxY == nExt[2])
1941 {
1942 ++maxY;
1943 }
1944 if (maxZ == nExt[4])
1945 {
1946 ++maxZ;
1947 }
1948 incY = extent[1] - extent[0];
1949 incZ = (extent[3] - extent[2]) * incY;
1950 outId = 0;
1951 inIdZ = incZ * (nExt[4] - extent[4]) + incY * (nExt[2] - extent[2]) + (nExt[0] - extent[0]);
1952 for (idxZ = nExt[4]; idxZ < maxZ; idxZ++)
1953 {
1954 inIdY = inIdZ;
1955 for (idxY = nExt[2]; idxY < maxY; idxY++)
1956 {
1957 inId = inIdY;
1958 for (idxX = nExt[0]; idxX < maxX; idxX++)
1959 {
1960 ncd->CopyData(this->CellData, inId, outId);
1961 ++inId;
1962 ++outId;
1963 }
1964 inIdY += incY;
1965 }
1966 inIdZ += incZ;
1967 }
1968
1969 this->PointData->ShallowCopy(npd);
1970 this->CellData->ShallowCopy(ncd);
1971 this->SetExtent(nExt);
1972 newImage->Delete();
1973 }
1974
1975 //------------------------------------------------------------------------------
GetScalarTypeMin(vtkInformation * meta_data)1976 double vtkImageData::GetScalarTypeMin(vtkInformation* meta_data)
1977 {
1978 return vtkDataArray::GetDataTypeMin(this->GetScalarType(meta_data));
1979 }
1980
1981 //------------------------------------------------------------------------------
GetScalarTypeMin()1982 double vtkImageData::GetScalarTypeMin()
1983 {
1984 return vtkDataArray::GetDataTypeMin(this->GetScalarType());
1985 }
1986
1987 //------------------------------------------------------------------------------
GetScalarTypeMax(vtkInformation * meta_data)1988 double vtkImageData::GetScalarTypeMax(vtkInformation* meta_data)
1989 {
1990 return vtkDataArray::GetDataTypeMax(this->GetScalarType(meta_data));
1991 }
1992
1993 //------------------------------------------------------------------------------
GetScalarTypeMax()1994 double vtkImageData::GetScalarTypeMax()
1995 {
1996 return vtkDataArray::GetDataTypeMax(this->GetScalarType());
1997 }
1998
1999 //------------------------------------------------------------------------------
SetExtent(int x1,int x2,int y1,int y2,int z1,int z2)2000 void vtkImageData::SetExtent(int x1, int x2, int y1, int y2, int z1, int z2)
2001 {
2002 int ext[6];
2003 ext[0] = x1;
2004 ext[1] = x2;
2005 ext[2] = y1;
2006 ext[3] = y2;
2007 ext[4] = z1;
2008 ext[5] = z2;
2009 this->SetExtent(ext);
2010 }
2011
2012 //------------------------------------------------------------------------------
SetDataDescription(int desc)2013 void vtkImageData::SetDataDescription(int desc)
2014 {
2015 if (desc == this->DataDescription)
2016 {
2017 return;
2018 }
2019
2020 this->DataDescription = desc;
2021
2022 if (this->Vertex)
2023 {
2024 this->Vertex->Delete();
2025 this->Vertex = nullptr;
2026 }
2027 if (this->Line)
2028 {
2029 this->Line->Delete();
2030 this->Line = nullptr;
2031 }
2032 if (this->Pixel)
2033 {
2034 this->Pixel->Delete();
2035 this->Pixel = nullptr;
2036 }
2037 if (this->Voxel)
2038 {
2039 this->Voxel->Delete();
2040 this->Voxel = nullptr;
2041 }
2042 switch (this->DataDescription)
2043 {
2044 case VTK_SINGLE_POINT:
2045 this->Vertex = vtkVertex::New();
2046 break;
2047
2048 case VTK_X_LINE:
2049 case VTK_Y_LINE:
2050 case VTK_Z_LINE:
2051 this->Line = vtkLine::New();
2052 break;
2053
2054 case VTK_XY_PLANE:
2055 case VTK_YZ_PLANE:
2056 case VTK_XZ_PLANE:
2057 this->Pixel = vtkPixel::New();
2058 break;
2059
2060 case VTK_XYZ_GRID:
2061 this->Voxel = vtkVoxel::New();
2062 break;
2063 }
2064 }
2065
2066 //------------------------------------------------------------------------------
SetExtent(int * extent)2067 void vtkImageData::SetExtent(int* extent)
2068 {
2069 int description;
2070
2071 description = vtkStructuredData::SetExtent(extent, this->Extent);
2072 if (description < 0) // improperly specified
2073 {
2074 vtkErrorMacro(<< "Bad Extent, retaining previous values");
2075 }
2076
2077 if (description == VTK_UNCHANGED)
2078 {
2079 return;
2080 }
2081
2082 vtkStructuredData::GetDimensionsFromExtent(extent, this->Dimensions);
2083
2084 this->SetDataDescription(description);
2085
2086 this->Modified();
2087 }
2088
2089 //------------------------------------------------------------------------------
GetDimensions()2090 int* vtkImageData::GetDimensions()
2091 {
2092 this->GetDimensions(this->Dimensions);
2093 return this->Dimensions;
2094 }
2095
2096 //------------------------------------------------------------------------------
GetDimensions(int * dOut)2097 void vtkImageData::GetDimensions(int* dOut)
2098 {
2099 const int* extent = this->Extent;
2100 dOut[0] = extent[1] - extent[0] + 1;
2101 dOut[1] = extent[3] - extent[2] + 1;
2102 dOut[2] = extent[5] - extent[4] + 1;
2103 }
2104
2105 #if VTK_ID_TYPE_IMPL != VTK_INT
2106 //------------------------------------------------------------------------------
GetDimensions(vtkIdType dims[3])2107 void vtkImageData::GetDimensions(vtkIdType dims[3])
2108 {
2109 // Use vtkIdType to avoid overflow on large images
2110 const int* extent = this->Extent;
2111 dims[0] = extent[1] - extent[0] + 1;
2112 dims[1] = extent[3] - extent[2] + 1;
2113 dims[2] = extent[5] - extent[4] + 1;
2114 }
2115 #endif
2116
2117 //------------------------------------------------------------------------------
SetAxisUpdateExtent(int idx,int min,int max,const int * updateExtent,int * axisUpdateExtent)2118 void vtkImageData::SetAxisUpdateExtent(
2119 int idx, int min, int max, const int* updateExtent, int* axisUpdateExtent)
2120 {
2121 if (idx > 2)
2122 {
2123 vtkWarningMacro("illegal axis!");
2124 return;
2125 }
2126
2127 memcpy(axisUpdateExtent, updateExtent, 6 * sizeof(int));
2128 if (axisUpdateExtent[idx * 2] != min)
2129 {
2130 axisUpdateExtent[idx * 2] = min;
2131 }
2132 if (axisUpdateExtent[idx * 2 + 1] != max)
2133 {
2134 axisUpdateExtent[idx * 2 + 1] = max;
2135 }
2136 }
2137
2138 //------------------------------------------------------------------------------
GetAxisUpdateExtent(int idx,int & min,int & max,const int * updateExtent)2139 void vtkImageData::GetAxisUpdateExtent(int idx, int& min, int& max, const int* updateExtent)
2140 {
2141 if (idx > 2)
2142 {
2143 vtkWarningMacro("illegal axis!");
2144 return;
2145 }
2146
2147 min = updateExtent[idx * 2];
2148 max = updateExtent[idx * 2 + 1];
2149 }
2150
2151 //------------------------------------------------------------------------------
GetActualMemorySize()2152 unsigned long vtkImageData::GetActualMemorySize()
2153 {
2154 return this->Superclass::GetActualMemorySize();
2155 }
2156
2157 //------------------------------------------------------------------------------
ShallowCopy(vtkDataObject * dataObject)2158 void vtkImageData::ShallowCopy(vtkDataObject* dataObject)
2159 {
2160 vtkImageData* imageData = vtkImageData::SafeDownCast(dataObject);
2161
2162 if (imageData != nullptr)
2163 {
2164 this->InternalImageDataCopy(imageData);
2165 }
2166
2167 // Do superclass
2168 this->Superclass::ShallowCopy(dataObject);
2169 }
2170
2171 //------------------------------------------------------------------------------
DeepCopy(vtkDataObject * dataObject)2172 void vtkImageData::DeepCopy(vtkDataObject* dataObject)
2173 {
2174 auto mkhold = vtkMemkindRAII(this->GetIsInMemkind());
2175 vtkImageData* imageData = vtkImageData::SafeDownCast(dataObject);
2176
2177 if (imageData != nullptr)
2178 {
2179 this->InternalImageDataCopy(imageData);
2180 }
2181
2182 // Do superclass
2183 this->Superclass::DeepCopy(dataObject);
2184 }
2185
2186 //------------------------------------------------------------------------------
2187 // This copies all the local variables (but not objects).
InternalImageDataCopy(vtkImageData * src)2188 void vtkImageData::InternalImageDataCopy(vtkImageData* src)
2189 {
2190 int idx;
2191
2192 // this->SetScalarType(src->GetScalarType());
2193 // this->SetNumberOfScalarComponents(src->GetNumberOfScalarComponents());
2194 for (idx = 0; idx < 3; ++idx)
2195 {
2196 this->Dimensions[idx] = src->Dimensions[idx];
2197 this->Increments[idx] = src->Increments[idx];
2198 this->Origin[idx] = src->Origin[idx];
2199 this->Spacing[idx] = src->Spacing[idx];
2200 }
2201 this->DirectionMatrix->DeepCopy(src->DirectionMatrix);
2202 this->ComputeTransforms();
2203 this->SetExtent(src->GetExtent());
2204 }
2205
2206 //------------------------------------------------------------------------------
GetNumberOfCells()2207 vtkIdType vtkImageData::GetNumberOfCells()
2208 {
2209 vtkIdType nCells = 1;
2210 int i;
2211 const int* extent = this->Extent;
2212
2213 vtkIdType dims[3];
2214 dims[0] = extent[1] - extent[0] + 1;
2215 dims[1] = extent[3] - extent[2] + 1;
2216 dims[2] = extent[5] - extent[4] + 1;
2217
2218 for (i = 0; i < 3; i++)
2219 {
2220 if (dims[i] == 0)
2221 {
2222 return 0;
2223 }
2224 if (dims[i] > 1)
2225 {
2226 nCells *= (dims[i] - 1);
2227 }
2228 }
2229
2230 return nCells;
2231 }
2232
2233 //============================================================================
2234 // Starting to make some more general methods that deal with any array
2235 // (not just scalars).
2236 //============================================================================
2237
2238 //------------------------------------------------------------------------------
2239 // This Method returns a pointer to a location in the vtkImageData.
2240 // Coordinates are in pixel units and are relative to the whole
2241 // image origin.
GetArrayIncrements(vtkDataArray * array,vtkIdType increments[3])2242 void vtkImageData::GetArrayIncrements(vtkDataArray* array, vtkIdType increments[3])
2243 {
2244 const int* extent = this->Extent;
2245 // We could store tuple increments and just
2246 // multiply by the number of components...
2247 increments[0] = array->GetNumberOfComponents();
2248 increments[1] = increments[0] * (extent[1] - extent[0] + 1);
2249 increments[2] = increments[1] * (extent[3] - extent[2] + 1);
2250 }
2251
2252 //------------------------------------------------------------------------------
GetArrayPointerForExtent(vtkDataArray * array,int extent[6])2253 void* vtkImageData::GetArrayPointerForExtent(vtkDataArray* array, int extent[6])
2254 {
2255 int tmp[3];
2256 tmp[0] = extent[0];
2257 tmp[1] = extent[2];
2258 tmp[2] = extent[4];
2259 return this->GetArrayPointer(array, tmp);
2260 }
2261
2262 //------------------------------------------------------------------------------
2263 // This Method returns am index to a location in the vtkImageData.
2264 // Coordinates are in pixel units and are relative to the whole
2265 // image origin.
GetTupleIndex(vtkDataArray * array,int coordinate[3])2266 vtkIdType vtkImageData::GetTupleIndex(vtkDataArray* array, int coordinate[3])
2267 {
2268 vtkIdType incs[3];
2269 vtkIdType idx;
2270
2271 if (array == nullptr)
2272 {
2273 return -1;
2274 }
2275
2276 const int* extent = this->Extent;
2277 // error checking: since most accesses will be from pointer arithmetic.
2278 // this should not waste much time.
2279 for (idx = 0; idx < 3; ++idx)
2280 {
2281 if (coordinate[idx] < extent[idx * 2] || coordinate[idx] > extent[idx * 2 + 1])
2282 {
2283 vtkErrorMacro(<< "GetPointer: Pixel (" << coordinate[0] << ", " << coordinate[1] << ", "
2284 << coordinate[2] << ") not in current extent: (" << extent[0] << ", "
2285 << extent[1] << ", " << extent[2] << ", " << extent[3] << ", " << extent[4]
2286 << ", " << extent[5] << ")");
2287 return -1;
2288 }
2289 }
2290
2291 // compute the index of the vector.
2292
2293 // Array increments incorporate the number of components, which is not how
2294 // vtkDataArrays are indexed. Instead, compute the tuple increments.
2295 {
2296 incs[0] = 1;
2297 incs[1] = (extent[1] - extent[0] + 1);
2298 incs[2] = incs[1] * (extent[3] - extent[2] + 1);
2299 }
2300
2301 idx = ((coordinate[0] - extent[0]) * incs[0] + (coordinate[1] - extent[2]) * incs[1] +
2302 (coordinate[2] - extent[4]) * incs[2]);
2303 // I could check to see if the array has the correct number
2304 // of tuples for the extent, but that would be an extra multiply.
2305 if (idx < 0 || idx > array->GetMaxId())
2306 {
2307 vtkErrorMacro("Coordinate (" << coordinate[0] << ", " << coordinate[1] << ", " << coordinate[2]
2308 << ") out side of array (max = " << array->GetMaxId());
2309 return -1;
2310 }
2311
2312 return idx;
2313 }
2314
2315 //------------------------------------------------------------------------------
2316 // This Method returns a pointer to a location in the vtkImageData.
2317 // Coordinates are in pixel units and are relative to the whole
2318 // image origin.
GetArrayPointer(vtkDataArray * array,int coordinate[3])2319 void* vtkImageData::GetArrayPointer(vtkDataArray* array, int coordinate[3])
2320 {
2321 return array->GetVoidPointer(
2322 array->GetNumberOfComponents() * this->GetTupleIndex(array, coordinate));
2323 }
2324
2325 //------------------------------------------------------------------------------
ComputeInternalExtent(int * intExt,int * tgtExt,int * bnds)2326 void vtkImageData::ComputeInternalExtent(int* intExt, int* tgtExt, int* bnds)
2327 {
2328 int i;
2329 const int* extent = this->Extent;
2330 for (i = 0; i < 3; ++i)
2331 {
2332 intExt[i * 2] = tgtExt[i * 2];
2333 if (intExt[i * 2] - bnds[i * 2] < extent[i * 2])
2334 {
2335 intExt[i * 2] = extent[i * 2] + bnds[i * 2];
2336 }
2337 intExt[i * 2 + 1] = tgtExt[i * 2 + 1];
2338 if (intExt[i * 2 + 1] + bnds[i * 2 + 1] > extent[i * 2 + 1])
2339 {
2340 intExt[i * 2 + 1] = extent[i * 2 + 1] - bnds[i * 2 + 1];
2341 }
2342 }
2343 }
2344
2345 //------------------------------------------------------------------------------
GetData(vtkInformation * info)2346 vtkImageData* vtkImageData::GetData(vtkInformation* info)
2347 {
2348 return info ? vtkImageData::SafeDownCast(info->Get(DATA_OBJECT())) : nullptr;
2349 }
2350
2351 //------------------------------------------------------------------------------
GetData(vtkInformationVector * v,int i)2352 vtkImageData* vtkImageData::GetData(vtkInformationVector* v, int i)
2353 {
2354 return vtkImageData::GetData(v->GetInformationObject(i));
2355 }
2356
2357 //------------------------------------------------------------------------------
SetSpacing(double i,double j,double k)2358 void vtkImageData::SetSpacing(double i, double j, double k)
2359 {
2360 vtkDebugMacro(<< this->GetClassName() << " (" << this << "): setting Spacing to (" << i << ","
2361 << j << "," << k << ")");
2362 if ((this->Spacing[0] != i) || (this->Spacing[1] != j) || (this->Spacing[2] != k))
2363 {
2364 this->Spacing[0] = i;
2365 this->Spacing[1] = j;
2366 this->Spacing[2] = k;
2367 this->ComputeTransforms();
2368 this->Modified();
2369 }
2370 }
2371
2372 //------------------------------------------------------------------------------
SetSpacing(const double ijk[3])2373 void vtkImageData::SetSpacing(const double ijk[3])
2374 {
2375 this->SetSpacing(ijk[0], ijk[1], ijk[2]);
2376 }
2377
2378 //------------------------------------------------------------------------------
SetOrigin(double i,double j,double k)2379 void vtkImageData::SetOrigin(double i, double j, double k)
2380 {
2381 vtkDebugMacro(<< this->GetClassName() << " (" << this << "): setting Origin to (" << i << "," << j
2382 << "," << k << ")");
2383 if ((this->Origin[0] != i) || (this->Origin[1] != j) || (this->Origin[2] != k))
2384 {
2385 this->Origin[0] = i;
2386 this->Origin[1] = j;
2387 this->Origin[2] = k;
2388 this->ComputeTransforms();
2389 this->Modified();
2390 }
2391 }
2392
2393 //------------------------------------------------------------------------------
SetOrigin(const double ijk[3])2394 void vtkImageData::SetOrigin(const double ijk[3])
2395 {
2396 this->SetOrigin(ijk[0], ijk[1], ijk[2]);
2397 }
2398
2399 //------------------------------------------------------------------------------
SetDirectionMatrix(vtkMatrix3x3 * m)2400 void vtkImageData::SetDirectionMatrix(vtkMatrix3x3* m)
2401 {
2402 vtkMTimeType lastModified = this->GetMTime();
2403 vtkSetObjectBodyMacro(DirectionMatrix, vtkMatrix3x3, m);
2404 if (lastModified < this->GetMTime())
2405 {
2406 this->ComputeTransforms();
2407 }
2408 }
2409
2410 //------------------------------------------------------------------------------
SetDirectionMatrix(const double elements[9])2411 void vtkImageData::SetDirectionMatrix(const double elements[9])
2412 {
2413 this->SetDirectionMatrix(elements[0], elements[1], elements[2], elements[3], elements[4],
2414 elements[5], elements[6], elements[7], elements[8]);
2415 }
2416
2417 //------------------------------------------------------------------------------
SetDirectionMatrix(double e00,double e01,double e02,double e10,double e11,double e12,double e20,double e21,double e22)2418 void vtkImageData::SetDirectionMatrix(double e00, double e01, double e02, double e10, double e11,
2419 double e12, double e20, double e21, double e22)
2420 {
2421 vtkMatrix3x3* m3 = this->DirectionMatrix;
2422 vtkMTimeType lastModified = m3->GetMTime();
2423
2424 m3->SetElement(0, 0, e00);
2425 m3->SetElement(0, 1, e01);
2426 m3->SetElement(0, 2, e02);
2427 m3->SetElement(1, 0, e10);
2428 m3->SetElement(1, 1, e11);
2429 m3->SetElement(1, 2, e12);
2430 m3->SetElement(2, 0, e20);
2431 m3->SetElement(2, 1, e21);
2432 m3->SetElement(2, 2, e22);
2433
2434 if (lastModified < m3->GetMTime())
2435 {
2436 this->ComputeTransforms();
2437 this->Modified();
2438 }
2439 }
2440
2441 //------------------------------------------------------------------------------
2442 template <typename T1, typename T2>
TransformCoordinates(T1 input0,T1 input1,T1 input2,T2 output[3],vtkMatrix4x4 * m4)2443 inline static void TransformCoordinates(
2444 T1 input0, T1 input1, T1 input2, T2 output[3], vtkMatrix4x4* m4)
2445 {
2446 double* mdata = m4->GetData();
2447 output[0] = mdata[0] * input0 + mdata[1] * input1 + mdata[2] * input2 + mdata[3];
2448 output[1] = mdata[4] * input0 + mdata[5] * input1 + mdata[6] * input2 + mdata[7];
2449 output[2] = mdata[8] * input0 + mdata[9] * input1 + mdata[10] * input2 + mdata[11];
2450 }
2451
2452 // must pass the inverse matrix
2453 template <typename T1, typename T2>
TransformNormal(T1 input0,T1 input1,T1 input2,T2 output[3],vtkMatrix4x4 * m4)2454 inline static void TransformNormal(T1 input0, T1 input1, T1 input2, T2 output[3], vtkMatrix4x4* m4)
2455 {
2456 double* mdata = m4->GetData();
2457 output[0] = mdata[0] * input0 + mdata[4] * input1 + mdata[8] * input2;
2458 output[1] = mdata[1] * input0 + mdata[5] * input1 + mdata[9] * input2;
2459 output[2] = mdata[2] * input0 + mdata[6] * input1 + mdata[10] * input2;
2460 }
2461
2462 // useful for when the ImageData is not available but the information
2463 // spacing, origin, direction are
TransformContinuousIndexToPhysicalPoint(double i,double j,double k,double const origin[3],double const spacing[3],double const direction[9],double xyz[3])2464 void vtkImageData::TransformContinuousIndexToPhysicalPoint(double i, double j, double k,
2465 double const origin[3], double const spacing[3], double const direction[9], double xyz[3])
2466 {
2467 for (int c = 0; c < 3; ++c)
2468 {
2469 xyz[c] = i * spacing[0] * direction[c * 3] + j * spacing[1] * direction[c * 3 + 1] +
2470 k * spacing[2] * direction[c * 3 + 2] + origin[c];
2471 }
2472 }
2473
2474 //------------------------------------------------------------------------------
TransformContinuousIndexToPhysicalPoint(double i,double j,double k,double xyz[3])2475 void vtkImageData::TransformContinuousIndexToPhysicalPoint(
2476 double i, double j, double k, double xyz[3])
2477 {
2478 TransformCoordinates<double, double>(i, j, k, xyz, this->IndexToPhysicalMatrix);
2479 }
2480
2481 //------------------------------------------------------------------------------
TransformContinuousIndexToPhysicalPoint(const double ijk[3],double xyz[3])2482 void vtkImageData::TransformContinuousIndexToPhysicalPoint(const double ijk[3], double xyz[3])
2483 {
2484
2485 TransformCoordinates<double, double>(ijk[0], ijk[1], ijk[2], xyz, this->IndexToPhysicalMatrix);
2486 }
2487
2488 //------------------------------------------------------------------------------
TransformIndexToPhysicalPoint(int i,int j,int k,double xyz[3])2489 void vtkImageData::TransformIndexToPhysicalPoint(int i, int j, int k, double xyz[3])
2490 {
2491 TransformCoordinates<int, double>(i, j, k, xyz, this->IndexToPhysicalMatrix);
2492 }
2493
2494 //------------------------------------------------------------------------------
TransformIndexToPhysicalPoint(const int ijk[3],double xyz[3])2495 void vtkImageData::TransformIndexToPhysicalPoint(const int ijk[3], double xyz[3])
2496 {
2497 TransformCoordinates<int, double>(ijk[0], ijk[1], ijk[2], xyz, this->IndexToPhysicalMatrix);
2498 }
2499
2500 //------------------------------------------------------------------------------
TransformPhysicalPointToContinuousIndex(double x,double y,double z,double ijk[3])2501 void vtkImageData::TransformPhysicalPointToContinuousIndex(
2502 double x, double y, double z, double ijk[3])
2503 {
2504 TransformCoordinates<double, double>(x, y, z, ijk, this->PhysicalToIndexMatrix);
2505 }
2506 //------------------------------------------------------------------------------
TransformPhysicalPointToContinuousIndex(const double xyz[3],double ijk[3])2507 void vtkImageData::TransformPhysicalPointToContinuousIndex(const double xyz[3], double ijk[3])
2508 {
2509 TransformCoordinates<double, double>(xyz[0], xyz[1], xyz[2], ijk, this->PhysicalToIndexMatrix);
2510 }
2511
2512 //------------------------------------------------------------------------------
TransformPhysicalNormalToContinuousIndex(const double xyz[3],double ijk[3])2513 void vtkImageData::TransformPhysicalNormalToContinuousIndex(const double xyz[3], double ijk[3])
2514 {
2515 TransformNormal<double, double>(xyz[0], xyz[1], xyz[2], ijk, this->IndexToPhysicalMatrix);
2516 }
2517
2518 //------------------------------------------------------------------------------
TransformPhysicalPlaneToContinuousIndex(double const normal[4],double xnormal[4])2519 void vtkImageData::TransformPhysicalPlaneToContinuousIndex(
2520 double const normal[4], double xnormal[4])
2521 {
2522 // transform the normal, note the inverse matrix is passed in
2523 TransformNormal<double, double>(
2524 normal[0], normal[1], normal[2], xnormal, this->IndexToPhysicalMatrix);
2525 vtkMath::Normalize(xnormal);
2526
2527 // transform the point
2528 double newPt[3];
2529 TransformCoordinates<double, double>(-normal[3] * normal[0], -normal[3] * normal[1],
2530 -normal[3] * normal[2], newPt, this->PhysicalToIndexMatrix);
2531
2532 // recompute plane eqn
2533 xnormal[3] = -xnormal[0] * newPt[0] - xnormal[1] * newPt[1] - xnormal[2] * newPt[2];
2534 }
2535
2536 //------------------------------------------------------------------------------
ComputeTransforms()2537 void vtkImageData::ComputeTransforms()
2538 {
2539 vtkMatrix4x4* m4 = vtkMatrix4x4::New();
2540 if (this->DirectionMatrix->IsIdentity())
2541 {
2542 m4->Zero();
2543 m4->SetElement(0, 0, this->Spacing[0]);
2544 m4->SetElement(1, 1, this->Spacing[1]);
2545 m4->SetElement(2, 2, this->Spacing[2]);
2546 m4->SetElement(3, 3, 1);
2547 }
2548 else
2549 {
2550 const double* m3 = this->DirectionMatrix->GetData();
2551 m4->SetElement(0, 0, m3[0] * this->Spacing[0]);
2552 m4->SetElement(0, 1, m3[1] * this->Spacing[1]);
2553 m4->SetElement(0, 2, m3[2] * this->Spacing[2]);
2554 m4->SetElement(1, 0, m3[3] * this->Spacing[0]);
2555 m4->SetElement(1, 1, m3[4] * this->Spacing[1]);
2556 m4->SetElement(1, 2, m3[5] * this->Spacing[2]);
2557 m4->SetElement(2, 0, m3[6] * this->Spacing[0]);
2558 m4->SetElement(2, 1, m3[7] * this->Spacing[1]);
2559 m4->SetElement(2, 2, m3[8] * this->Spacing[2]);
2560 m4->SetElement(3, 0, 0);
2561 m4->SetElement(3, 1, 0);
2562 m4->SetElement(3, 2, 0);
2563 m4->SetElement(3, 3, 1);
2564 }
2565 m4->SetElement(0, 3, this->Origin[0]);
2566 m4->SetElement(1, 3, this->Origin[1]);
2567 m4->SetElement(2, 3, this->Origin[2]);
2568
2569 this->IndexToPhysicalMatrix->DeepCopy(m4);
2570 vtkMatrix4x4::Invert(m4, this->PhysicalToIndexMatrix);
2571 m4->Delete();
2572 }
2573
2574 //------------------------------------------------------------------------------
ComputeIndexToPhysicalMatrix(double const origin[3],double const spacing[3],double const direction[9],double result[16])2575 void vtkImageData::ComputeIndexToPhysicalMatrix(
2576 double const origin[3], double const spacing[3], double const direction[9], double result[16])
2577 {
2578 for (int i = 0; i < 3; ++i)
2579 {
2580 result[i * 4] = direction[i * 3] * spacing[0];
2581 result[i * 4 + 1] = direction[i * 3 + 1] * spacing[1];
2582 result[i * 4 + 2] = direction[i * 3 + 2] * spacing[2];
2583 }
2584
2585 result[3] = origin[0];
2586 result[7] = origin[1];
2587 result[11] = origin[2];
2588 result[12] = 0;
2589 result[13] = 0;
2590 result[14] = 0;
2591 result[15] = 1;
2592 }
2593
2594 //------------------------------------------------------------------------------
HasAnyBlankPoints()2595 bool vtkImageData::HasAnyBlankPoints()
2596 {
2597 return this->IsAnyBitSet(this->GetPointGhostArray(), vtkDataSetAttributes::HIDDENPOINT);
2598 }
2599
2600 //------------------------------------------------------------------------------
HasAnyBlankCells()2601 bool vtkImageData::HasAnyBlankCells()
2602 {
2603 int cellBlanking = this->IsAnyBitSet(this->GetCellGhostArray(), vtkDataSetAttributes::HIDDENCELL);
2604 return cellBlanking || this->HasAnyBlankPoints();
2605 }
2606