Updates will be applied today, at 4 pm EDT (UTC-0400). No interruption, and site should remain available.

vtkVolumeOutlineSource.cxx 22.5 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkVolumeOutlineSource.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/
#include "vtkVolumeOutlineSource.h"

17 18 19
#include "vtkCellArray.h"
#include "vtkCellArrayIterator.h"
#include "vtkCellData.h"
20
#include "vtkDataSet.h"
21
#include "vtkDemandDrivenPipeline.h"
22
#include "vtkIdList.h"
23
#include "vtkImageData.h"
24 25
#include "vtkInformation.h"
#include "vtkInformationVector.h"
26
#include "vtkMath.h"
27
#include "vtkNew.h"
28
#include "vtkObjectFactory.h"
29
#include "vtkPoints.h"
30
#include "vtkPolyData.h"
31
#include "vtkStreamingDemandDrivenPipeline.h"
32 33 34 35 36
#include "vtkUnsignedCharArray.h"
#include "vtkVolumeMapper.h"

vtkStandardNewMacro(vtkVolumeOutlineSource);

37
vtkCxxSetObjectMacro(vtkVolumeOutlineSource, VolumeMapper, vtkVolumeMapper);
38 39

//----------------------------------------------------------------------------
40
vtkVolumeOutlineSource::vtkVolumeOutlineSource()
41
{
42
  this->VolumeMapper = nullptr;
43
  this->GenerateScalars = 0;
44
  this->GenerateOutline = 1;
45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
  this->GenerateFaces = 0;
  this->ActivePlaneId = -1;

  this->Color[0] = 1.0;
  this->Color[1] = 0.0;
  this->Color[2] = 0.0;

  this->ActivePlaneColor[0] = 1.0;
  this->ActivePlaneColor[1] = 1.0;
  this->ActivePlaneColor[2] = 0.0;

  this->SetNumberOfInputPorts(0);
}

//----------------------------------------------------------------------------
60
vtkVolumeOutlineSource::~vtkVolumeOutlineSource()
61 62
{
  if (this->VolumeMapper)
63
  {
64
    this->VolumeMapper->Delete();
65
    this->VolumeMapper = nullptr;
66
  }
67 68 69 70 71
}

//----------------------------------------------------------------------------
void vtkVolumeOutlineSource::PrintSelf(ostream& os, vtkIndent indent)
{
72
  this->Superclass::PrintSelf(os, indent);
73 74 75

  os << indent << "VolumeMapper: ";
  if (this->VolumeMapper)
76
  {
77
    os << this->VolumeMapper << "\n";
78
  }
79
  else
80
  {
81
    os << "(none)\n";
82
  }
83

84
  os << indent << "GenerateFaces: " << (this->GenerateFaces ? "On\n" : "Off\n");
85

86
  os << indent << "GenerateOutline: " << (this->GenerateOutline ? "On\n" : "Off\n");
87

88
  os << indent << "GenerateScalars: " << (this->GenerateScalars ? "On\n" : "Off\n");
89

90 91
  os << indent << "Color: " << this->Color[0] << ", " << this->Color[1] << ", " << this->Color[2]
     << "\n";
92 93 94 95 96 97 98 99 100 101 102

  os << indent << "ActivePlaneId: " << this->ActivePlaneId << "\n";

  os << indent << "ActivePlaneColor: " << this->ActivePlaneColor[0] << ", "
     << this->ActivePlaneColor[1] << ", " << this->ActivePlaneColor[2] << "\n";
}

//----------------------------------------------------------------------------
int vtkVolumeOutlineSource::ComputeCubePlanes(
  double planes[3][4], double croppingPlanes[6], double bounds[6])
{
103 104 105 106 107 108
  // Combine the CroppingRegionPlanes and the Bounds to create
  // a single array.  For each dimension, store the planes in
  // the following order: lo_bound, lo_crop_plane, hi_crop_plane, hi_bound.
  // Also do range checking to ensure that the cropping planes
  // are clamped to the bound limits.

109
  for (int i = 0; i < 3; i++)
110
  {
111 112
    int j0 = 2 * i;
    int j1 = 2 * i + 1;
113 114 115 116 117 118

    double a = bounds[j0];
    double b = croppingPlanes[j0];
    double c = croppingPlanes[j1];
    double d = bounds[j1];

119
    // Sanity check
120
    if (a > d || b > c)
121
    {
122
      return 0;
123
    }
124

125
    // Clamp cropping planes to bounds
126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141
    if (b < a)
    {
      b = a;
    }
    if (b > d)
    {
      b = d;
    }
    if (c < a)
    {
      c = a;
    }
    if (c > d)
    {
      c = d;
    }
142 143 144 145 146

    planes[i][0] = a;
    planes[i][1] = b;
    planes[i][2] = c;
    planes[i][3] = d;
147
  }
148 149 150 151 152

  return 1;
}

//----------------------------------------------------------------------------
153 154 155
int vtkVolumeOutlineSource::ComputePipelineMTime(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** vtkNotUsed(inputVector), vtkInformationVector* vtkNotUsed(outputVector),
  int vtkNotUsed(requestFromOutputPort), vtkMTimeType* mtime)
156
{
Bill Lorensen's avatar
Bill Lorensen committed
157
  vtkMTimeType mTime = this->GetMTime();
158
  if (this->VolumeMapper)
159
  {
Bill Lorensen's avatar
Bill Lorensen committed
160
    vtkMTimeType mapperMTime = this->VolumeMapper->GetMTime();
161
    if (mapperMTime > mTime)
162
    {
163
      mTime = mapperMTime;
164
    }
165 166
    vtkDemandDrivenPipeline* input =
      vtkDemandDrivenPipeline::SafeDownCast(this->VolumeMapper->GetInputExecutive());
167
    if (input)
168
    {
169 170 171
      // Need to do this because we are not formally connected
      // to the Mapper's pipeline
      input->UpdateInformation();
Bill Lorensen's avatar
Bill Lorensen committed
172
      vtkMTimeType pipelineMTime = input->GetPipelineMTime();
173
      if (pipelineMTime > mTime)
174
      {
175 176 177
        mTime = pipelineMTime;
      }
    }
178
  }
179 180 181 182 183 184 185

  *mtime = mTime;

  return 1;
}

//----------------------------------------------------------------------------
186 187
int vtkVolumeOutlineSource::RequestInformation(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** vtkNotUsed(inputVector), vtkInformationVector* vtkNotUsed(outputVector))
188 189 190 191 192
{
  // Get the mapper's input, since this is the most convenient
  // place to do so.

  if (!this->VolumeMapper)
193
  {
194 195
    vtkWarningMacro("No VolumeMapper has been set.");
    return 1;
196
  }
197

198
  vtkInformation* mapInfo = this->VolumeMapper->GetInputInformation();
199

200
  if (!mapInfo)
201
  {
202 203
    vtkWarningMacro("The VolumeMapper does not have an input set.");
    return 1;
204
  }
205 206 207 208 209

  // Don't have to update mapper's input, since it was done in
  // ComputePipelineMTime.
  // data->UpdateInformation();

210 211
  // Don't call GetBounds because we need WholeExtent, while
  // GetBounds only returns the bounds for Extent.
212 213 214 215 216

  double spacing[3];
  double origin[3];
  int extent[6];

217 218
  mapInfo->Get(vtkDataObject::SPACING(), spacing);
  mapInfo->Get(vtkDataObject::ORIGIN(), origin);
219
  mapInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), extent);
220 221

  for (int i = 0; i < 3; i++)
222
  {
223 224
    int j0 = 2 * i;
    int j1 = j0 + 1;
225 226

    if (extent[j0] > extent[j1])
227
    {
228 229
      vtkMath::UninitializeBounds(this->Bounds);
      break;
230
    }
231 232

    if (spacing[i] > 0)
233
    {
234 235
      this->Bounds[j0] = origin[i] + spacing[i] * extent[j0];
      this->Bounds[j1] = origin[i] + spacing[i] * extent[j1];
236
    }
237
    else
238
    {
239 240
      this->Bounds[j0] = origin[i] + spacing[i] * extent[j1];
      this->Bounds[j1] = origin[i] + spacing[i] * extent[j0];
241
    }
242 243 244

    this->CroppingRegionPlanes[j0] = this->Bounds[j0];
    this->CroppingRegionPlanes[j1] = this->Bounds[j1];
245
  }
246 247 248 249 250

  this->CroppingRegionFlags = 0x0002000;

  this->Cropping = this->VolumeMapper->GetCropping();
  if (this->Cropping)
251
  {
252 253
    this->CroppingRegionFlags = this->VolumeMapper->GetCroppingRegionFlags();
    this->VolumeMapper->GetCroppingRegionPlanes(this->CroppingRegionPlanes);
254
  }
255 256 257 258 259

  return 1;
}

//----------------------------------------------------------------------------
260 261
int vtkVolumeOutlineSource::RequestData(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** vtkNotUsed(inputVector), vtkInformationVector* outputVector)
262 263
{
  // get the info object
264
  vtkInformation* outInfo = outputVector->GetInformationObject(0);
265 266

  // get the output
267
  vtkPolyData* output = vtkPolyData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
268 269 270

  vtkDebugMacro(<< "Creating cropping region outline");

271 272
  // For each of the 3 dimensions, there are 4 planes: two bounding planes
  // on the outside, and two cropping region planes inside.
273 274 275
  double planes[3][4];

  if (!this->VolumeMapper || !this->VolumeMapper->GetInput() ||
276
    !this->ComputeCubePlanes(planes, this->CroppingRegionPlanes, this->Bounds))
277
  {
278
    // If the bounds or the cropping planes are invalid, clear the data
279 280 281
    output->SetPoints(nullptr);
    output->SetLines(nullptr);
    output->GetCellData()->SetScalars(nullptr);
282 283

    return 1;
284
  }
285

286
  // Compute the tolerance for considering points or planes to be coincident
287 288
  double tol = 0;
  for (int planeDim = 0; planeDim < 3; planeDim++)
289
  {
290
    double d = planes[planeDim][3] - planes[planeDim][0];
291
    tol += d * d;
292
  }
293
  tol = sqrt(tol) * 1e-5;
294

295 296
  // Create an array to nudge crop planes over to the bounds if they are
  // within tolerance of the bounds
297 298 299 300 301 302
  int tolPtId[3][4];
  this->NudgeCropPlanesToBounds(tolPtId, planes, tol);

  // The all-important cropping flags
  int flags = this->CroppingRegionFlags;

303
  // The active plane, which gets a special color for its scalars
304
  int activePlane = this->ActivePlaneId;
305 306 307 308
  if (activePlane > 5)
  {
    activePlane = -1;
  }
309 310 311 312 313 314

  // Convert the colors to unsigned char for scalars
  unsigned char colors[2][3];
  this->CreateColorValues(colors, this->Color, this->ActivePlaneColor);

  // Create the scalars used to color the lines
315
  vtkUnsignedCharArray* scalars = nullptr;
316 317

  if (this->GenerateScalars)
318
  {
319 320
    scalars = vtkUnsignedCharArray::New();
    scalars->SetNumberOfComponents(3);
321
  }
322 323

  // Generate all the lines for the outline.
324
  vtkCellArray* lines = nullptr;
325 326

  if (this->GenerateOutline)
327
  {
328 329
    lines = vtkCellArray::New();
    this->GenerateLines(lines, scalars, colors, activePlane, flags, tolPtId);
330
  }
331 332

  // Generate the polys for the outline
333
  vtkCellArray* polys = nullptr;
334 335

  if (this->GenerateFaces)
336
  {
337 338
    polys = vtkCellArray::New();
    this->GeneratePolys(polys, scalars, colors, activePlane, flags, tolPtId);
339
  }
340 341

  // Generate the points that are used by the lines.
342
  vtkPoints* points = vtkPoints::New();
343 344 345 346 347 348 349
  this->GeneratePoints(points, lines, polys, planes, tol);

  output->SetPoints(points);
  points->Delete();

  output->SetPolys(polys);
  if (polys)
350
  {
351
    polys->Delete();
352
  }
353 354

  output->SetLines(lines);
355
  if (lines)
356
  {
357
    lines->Delete();
358
  }
359 360 361

  output->GetCellData()->SetScalars(scalars);
  if (scalars)
362
  {
363
    scalars->Delete();
364
  }
365 366 367 368 369

  return 1;
}

//----------------------------------------------------------------------------
370 371
void vtkVolumeOutlineSource::GeneratePolys(vtkCellArray* polys, vtkUnsignedCharArray* scalars,
  unsigned char colors[2][3], int activePlane, int flags, int tolPtId[3][4])
372
{
373
  // Loop over the three dimensions and create the face rectangles
374
  for (int dim0 = 0; dim0 < 3; dim0++)
375
  {
376
    // Compute the other two dimension indices
377 378
    int dim1 = (dim0 + 1) % 3;
    int dim2 = (dim0 + 2) % 3;
379 380 381 382 383 384

    // Indices into the cubes
    int idx[3];

    // Loop over the "dim+2" dimension
    for (int i = 0; i < 4; i++)
385
    {
386 387 388 389
      idx[dim2] = i;

      // Loop over the "dim+1" dimension
      for (int j = 0; j < 3; j++)
390
      {
391 392 393
        idx[dim1] = j;

        // Make sure that the rect dim is not less than tolerance
394
        if ((j == 0 && tolPtId[dim1][1] == 0) || (j == 2 && tolPtId[dim1][2] == 3))
395
        {
396
          continue;
397
        }
398 399 400

        // Loop over rectangle along the "dim" dimension
        for (int k = 0; k < 3; k++)
401
        {
402 403 404
          idx[dim0] = k;

          // Make sure that the rect dim is not less than tolerance
405
          if ((k == 0 && tolPtId[dim0][1] == 0) || (k == 2 && tolPtId[dim0][2] == 3))
406
          {
407
            continue;
408
          }
409 410 411 412 413

          // The points in the rectangle, which are nudged over to the
          // volume bounds if the cropping planes are within tolerance
          // of the volume bounds.
          int pointId[4];
414
          pointId[0] = (tolPtId[2][idx[2]] * 16 + tolPtId[1][idx[1]] * 4 + tolPtId[0][idx[0]]);
415
          idx[dim0] = k + 1;
416
          pointId[1] = (tolPtId[2][idx[2]] * 16 + tolPtId[1][idx[1]] * 4 + tolPtId[0][idx[0]]);
417
          idx[dim1] = j + 1;
418
          pointId[2] = (tolPtId[2][idx[2]] * 16 + tolPtId[1][idx[1]] * 4 + tolPtId[0][idx[0]]);
419
          idx[dim0] = k;
420
          pointId[3] = (tolPtId[2][idx[2]] * 16 + tolPtId[1][idx[1]] * 4 + tolPtId[0][idx[0]]);
421 422 423 424
          idx[dim1] = j;

          // Loop through the two cubes adjacent to the rectangle,
          // in order to determine whether the rectangle is internal:
425 426
          // only external faces will be drawn.  The "bitCheck"
          // holds a bit for each of these two cubes.
427 428 429 430 431
          int bitCheck = 0;
          int cidx[3];
          cidx[dim0] = idx[dim0];
          cidx[dim1] = idx[dim1];
          for (int ii = 0; ii < 2; ii++)
432
          {
433 434 435 436
            // First get idx[dim2]-1, then idx[dim2]
            cidx[dim2] = idx[dim2] + ii - 1;
            int flagval = 0;
            if (cidx[dim2] >= 0 && cidx[dim2] < 3)
437
            {
438
              int flagbit = cidx[2] * 9 + cidx[1] * 3 + cidx[0];
439
              flagval = ((flags >> flagbit) & 1);
440
            }
441 442
            bitCheck <<= 1;
            bitCheck |= flagval;
443
          }
444

445 446
          // Whether we need to create a face depends on bitCheck.
          // Values 00, 11 don't need lines, while 01 and 10 do.
447 448 449

          // If our rect isn't an internal rect
          if (bitCheck != 0x0 && bitCheck != 0x3)
450
          {
451 452 453
            // Check if the rect is on our active plane
            int active = 0;
            if (activePlane >= 0)
454
            {
455
              int planeDim = (activePlane >> 1);    // same as "/ 2"
456 457
              int planeIdx = 1 + (activePlane & 1); // same as "% 2"
              if (planeDim == dim2 && i == planeIdx)
458
              {
459 460
                active = 1;
              }
461
            }
462 463 464 465

            // Insert the rectangle with the correct sense
            polys->InsertNextCell(4);
            if (bitCheck == 0x2)
466
            {
467 468 469 470
              polys->InsertCellPoint(pointId[0]);
              polys->InsertCellPoint(pointId[1]);
              polys->InsertCellPoint(pointId[2]);
              polys->InsertCellPoint(pointId[3]);
471
            }
472
            else // (bitCheck == 0x1)
473
            {
474 475 476 477
              polys->InsertCellPoint(pointId[3]);
              polys->InsertCellPoint(pointId[2]);
              polys->InsertCellPoint(pointId[1]);
              polys->InsertCellPoint(pointId[0]);
478
            }
479

480
            // Color the face
481
            if (scalars)
482
            {
David C. Lonie's avatar
David C. Lonie committed
483
              scalars->InsertNextTypedTuple(colors[active]);
484
            }
485
          }
486

487
        } // loop over k
488 489 490
      }   // loop over j
    }     // loop over i
  }       // loop over dim0
491 492 493
}

//----------------------------------------------------------------------------
494 495
void vtkVolumeOutlineSource::GenerateLines(vtkCellArray* lines, vtkUnsignedCharArray* scalars,
  unsigned char colors[2][3], int activePlane, int flags, int tolPtId[3][4])
496 497 498
{
  // Loop over the three dimensions and create the lines
  for (int dim0 = 0; dim0 < 3; dim0++)
499
  {
500
    // Compute the other two dimension indices
501 502
    int dim1 = (dim0 + 1) % 3;
    int dim2 = (dim0 + 2) % 3;
503 504 505 506 507 508

    // Indices into the cubes
    int idx[3];

    // Loop over the "dim+2" dimension
    for (int i = 0; i < 4; i++)
509
    {
510 511 512 513
      idx[dim2] = i;

      // Loop over the "dim+1" dimension
      for (int j = 0; j < 4; j++)
514
      {
515 516 517 518
        idx[dim1] = j;

        // Loop over line segments along the "dim" dimension
        for (int k = 0; k < 3; k++)
519
        {
520 521 522
          idx[dim0] = k;

          // Make sure that the segment length is not less than tolerance
523
          if ((k == 0 && tolPtId[dim0][1] == 0) || (k == 2 && tolPtId[dim0][2] == 3))
524
          {
525
            continue;
526
          }
527 528 529 530

          // The endpoints of the segment, which are nudged over to the
          // volume bounds if the cropping planes are within tolerance
          // of the volume bounds.
531
          int pointId0 = (tolPtId[2][idx[2]] * 16 + tolPtId[1][idx[1]] * 4 + tolPtId[0][idx[0]]);
532
          idx[dim0] = k + 1;
533
          int pointId1 = (tolPtId[2][idx[2]] * 16 + tolPtId[1][idx[1]] * 4 + tolPtId[0][idx[0]]);
534 535 536
          idx[dim0] = k;

          // Loop through the four cubes adjacent to the line segment,
537
          // in order to determine whether the line segment is on an
538 539 540 541 542 543
          // edge: only the edge lines will be drawn.  The "bitCheck"
          // holds a bit for each of these four cubes.
          int bitCheck = 0;
          int cidx[3];
          cidx[dim0] = idx[dim0];
          for (int ii = 0; ii < 2; ii++)
544
          {
545 546 547
            // First get idx[dim1]-1, then idx[dim1]
            cidx[dim1] = idx[dim1] + ii - 1;
            for (int jj = 0; jj < 2; jj++)
548
            {
549 550
              // First get idx[dim2]-1, then idx[dim2], but reverse
              // the order when ii loop is on its second iteration
551
              cidx[dim2] = idx[dim2] + (ii ^ jj) - 1;
552
              int flagval = 0;
553
              if (cidx[dim1] >= 0 && cidx[dim1] < 3 && cidx[dim2] >= 0 && cidx[dim2] < 3)
554
              {
555
                int flagbit = cidx[2] * 9 + cidx[1] * 3 + cidx[0];
556
                flagval = ((flags >> flagbit) & 1);
557
              }
558 559 560
              bitCheck <<= 1;
              bitCheck |= flagval;
            }
561
          }
562

luz.paz's avatar
luz.paz committed
563
          // Whether we need a line depends on the value of bitCheck.
564
          // Values 0000, 0011, 0110, 1100, 1001, 1111 don't need lines.
565
          // Build a bitfield to check our bitfield values against, each
566
          // set bit in this new bitfield corresponds to a non-edge case.
567 568
          const int noLineValues =
            ((1 << 0x0) | (1 << 0x3) | (1 << 0x6) | (1 << 0x9) | (1 << 0xc) | (1 << 0xf));
569 570 571

          // If our line segment is an edge, there is lots of work to do.
          if (((noLineValues >> bitCheck) & 1) == 0)
572
          {
573 574 575
            // Check if the line segment is on our active plane
            int active = 0;
            if (activePlane >= 0)
576
            {
577
              int planeDim = (activePlane >> 1);    // same as "/ 2"
578
              int planeIdx = 1 + (activePlane & 1); // same as "% 2"
579
              if ((planeDim == dim2 && i == planeIdx) || (planeDim == dim1 && j == planeIdx))
580
              {
581 582
                active = 1;
              }
583
            }
584 585 586 587

            // Check to make sure line segment isn't already there
            int foundDuplicate = 0;
            lines->InitTraversal();
588
            vtkIdType npts;
589
            const vtkIdType* pts;
590
            for (int cellId = 0; lines->GetNextCell(npts, pts); cellId++)
591
            {
592
              if (pts[0] == pointId0 && pts[1] == pointId1)
593
              {
594 595
                // Change color if current segment is on active plane
                if (scalars && active)
596
                {
David C. Lonie's avatar
David C. Lonie committed
597
                  scalars->SetTypedTuple(cellId, colors[active]);
598
                }
599 600 601
                foundDuplicate = 1;
                break;
              }
602
            }
603 604

            if (!foundDuplicate)
605
            {
606 607 608 609 610 611 612
              // Insert the line segment
              lines->InsertNextCell(2);
              lines->InsertCellPoint(pointId0);
              lines->InsertCellPoint(pointId1);

              // Color the line segment
              if (scalars)
613
              {
David C. Lonie's avatar
David C. Lonie committed
614
                scalars->InsertNextTypedTuple(colors[active]);
615
              }
616
            }
617
          }
618

619
        } // loop over k
620 621 622
      }   // loop over j
    }     // loop over i
  }       // loop over dim0
623 624 625 626
}

//----------------------------------------------------------------------------
void vtkVolumeOutlineSource::GeneratePoints(
627
  vtkPoints* points, vtkCellArray* lines, vtkCellArray* polys, double planes[3][4], double tol)
628 629 630 631 632 633
{
  // Use a bitfield to store which of the 64 points we need.
  // Two 32-bit ints are a convenient, portable way to do this.
  unsigned int pointBits1 = 0;
  unsigned int pointBits2 = 0;

634
  vtkIdType npts;
635 636
  const vtkIdType* pts;
  vtkCellArray* cellArrays[2];
637 638
  cellArrays[0] = lines;
  cellArrays[1] = polys;
639

640
  for (int arrayId = 0; arrayId < 2; arrayId++)
641
  {
642
    if (cellArrays[arrayId])
643
    {
644 645
      cellArrays[arrayId]->InitTraversal();
      while (cellArrays[arrayId]->GetNextCell(npts, pts))
646
      {
647
        for (int ii = 0; ii < npts; ii++)
648
        {
649
          int pointId = pts[ii];
650 651 652 653 654 655 656 657
          if (pointId < 32)
          {
            pointBits1 |= (1 << pointId);
          }
          else
          {
            pointBits2 |= (1 << (pointId - 32));
          }
658 659 660
        }
      }
    }
661
  }
662 663 664 665 666 667 668 669

  // Create the array of up to 64 points, and use the pointBits bitfield
  // to find out which points were used.  It is also necessary to go through
  // and update the cells with the modified point ids.
  unsigned int pointBits = pointBits1;
  int ptId = 0;
  int newPtId = 0;

670
  vtkNew<vtkIdList> repCell;
671
  for (int i = 0; i < 4; i++)
672
  {
673
    // If we're halfway done, switch over to the next 32 bits
674 675 676 677
    if (i == 2)
    {
      pointBits = pointBits2;
    }
678

679
    for (int j = 0; j < 4; j++)
680
    {
681
      for (int k = 0; k < 4; k++)
682
      {
683
        // Check to see if this point was actually used
684
        if ((pointBits & 1))
685
        {
686
          // Add or subtract tolerance as an offset to help depth check
687 688 689
          double x = planes[0][k] + tol * (1 - 2 * (k < 2));
          double y = planes[1][j] + tol * (1 - 2 * (j < 2));
          double z = planes[2][i] + tol * (1 - 2 * (i < 2));
690 691 692

          points->InsertNextPoint(x, y, z);

693
          for (int arrayId = 0; arrayId < 2; arrayId++)
694
          {
695
            // Go through the cells, substitute old Id for new Id
696
            vtkCellArray* cells = cellArrays[arrayId];
697
            if (cells)
698
            {
699
              auto cellIter = vtk::TakeSmartPointer(cells->NewIterator());
700
              for (cellIter->GoToFirstCell(); !cellIter->IsDoneWithTraversal();
701
                   cellIter->GoToNextCell())
702
              {
703 704
                cellIter->GetCurrentCell(repCell);
                for (int ii = 0; ii < repCell->GetNumberOfIds(); ii++)
705
                {
706 707 708 709
                  if (repCell->GetId(ii) == ptId)
                  {
                    repCell->SetId(ii, newPtId);
                  }
710
                }
711
                cellIter->ReplaceCurrentCell(repCell);
712 713 714
              }
            }
          }
715 716
          newPtId++;
        }
717 718 719 720
        pointBits >>= 1;
        ptId++;
      }
    }
721
  }
722 723 724 725 726 727 728
}

//----------------------------------------------------------------------------
void vtkVolumeOutlineSource::NudgeCropPlanesToBounds(
  int tolPtId[3][4], double planes[3][4], double tol)
{
  for (int dim = 0; dim < 3; dim++)
729
  {
730 731 732 733 734 735 736 737 738 739 740 741
    tolPtId[dim][0] = 0;
    tolPtId[dim][1] = 1;
    tolPtId[dim][2] = 2;
    tolPtId[dim][3] = 3;
    if (planes[dim][1] - planes[dim][0] < tol)
    {
      tolPtId[dim][1] = 0;
    }
    if (planes[dim][3] - planes[dim][2] < tol)
    {
      tolPtId[dim][2] = 3;
    }
742
  }
743 744 745 746 747 748 749
}

//----------------------------------------------------------------------------
void vtkVolumeOutlineSource::CreateColorValues(
  unsigned char colors[2][3], double color1[3], double color2[3])
{
  // Convert the two colors to unsigned char
750
  double* dcolors[2];
751 752 753 754
  dcolors[0] = color1;
  dcolors[1] = color2;

  for (int i = 0; i < 2; i++)
755
  {
756
    for (int j = 0; j < 3; j++)
757
    {
758
      double val = dcolors[i][j];
759 760 761 762 763 764 765 766 767
      if (val < 0)
      {
        val = 0;
      }
      if (val > 1)
      {
        val = 1;
      }
      colors[i][j] = static_cast<unsigned char>(val * 255);
768
    }
769
  }
770
}