vtkSlicer.C 25.4 KB
Newer Older
hrchilds's avatar
hrchilds committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    $RCSfile: vtkSlicer.cxx,v $
  Language:  C++
  Date:      $Date: 2002/02/22 21:16:54 $
  Version:   $Revision: 1.66 $

  Copyright (c) 1993-2002 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 "vtkSlicer.h"
hrchilds's avatar
hrchilds committed
19

20 21 22
#include <math.h>
#include <vector>

hrchilds's avatar
hrchilds committed
23
#include <vtkAppendPolyData.h>
hrchilds's avatar
hrchilds committed
24
#include <vtkCellData.h>
25
#include <vtkCutter.h>
hrchilds's avatar
hrchilds committed
26
#include <vtkDataSet.h>
hrchilds's avatar
hrchilds committed
27
#include <vtkFloatArray.h>
28 29
#include <vtkInformation.h>
#include <vtkInformationVector.h>
hrchilds's avatar
hrchilds committed
30
#include <vtkObjectFactory.h>
hrchilds's avatar
hrchilds committed
31 32 33
#include <vtkPlane.h>
#include <vtkPointData.h>
#include <vtkPolyData.h>
hrchilds's avatar
hrchilds committed
34 35
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>
hrchilds's avatar
hrchilds committed
36 37
#include <vtkSurfaceFromVolume.h>
#include <vtkTriangulationTables.h>
hrchilds's avatar
hrchilds committed
38 39
#include <vtkUnstructuredGrid.h>

40
#include <vtkCreateTriangleHelpers.h>
41 42
#include <vtkVisItCutter.h>

hrchilds's avatar
hrchilds committed
43 44 45 46 47 48 49 50 51 52 53 54 55 56
vtkStandardNewMacro(vtkSlicer);

vtkSlicer::vtkSlicer()
{
  this->CellList = NULL;
  this->CellListSize = 0;
  Origin[0] = Origin[1] = Origin[2] = 0.;
  Normal[0] = Normal[1] = Normal[2] = 0.;
}

vtkSlicer::~vtkSlicer()
{
}

57 58
void
vtkSlicer::SetCellList(vtkIdType *cl, vtkIdType size)
hrchilds's avatar
hrchilds committed
59 60 61 62 63
{
    this->CellList = cl;
    this->CellListSize = size;
}

64 65 66 67 68 69 70 71 72 73 74 75 76 77
// ****************************************************************************
//  Method: vtkSlicer::RequestData.
//
//  Modifications:
//    Eric Brugger, Thu Jan 10 10:24:20 PST 2013
//    Modified to inherit from vtkPolyDataAlgorithm.
//
// ****************************************************************************

int
vtkSlicer::RequestData(
    vtkInformation *vtkNotUsed(request),
    vtkInformationVector **inputVector,
    vtkInformationVector *outputVector)
hrchilds's avatar
hrchilds committed
78
{
79 80 81 82 83 84 85 86 87 88 89 90 91
    vtkDebugMacro(<<"Executing vtkSlicer");

    // get the info objects
    vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
    vtkInformation *outInfo = outputVector->GetInformationObject(0);

    //
    // Initialize some frequently used values.
    //
    input  = vtkDataSet::SafeDownCast(
        inInfo->Get(vtkDataObject::DATA_OBJECT()));
    output = vtkPolyData::SafeDownCast(
        outInfo->Get(vtkDataObject::DATA_OBJECT()));
hrchilds's avatar
hrchilds committed
92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109

    int do_type = input->GetDataObjectType();
    if (do_type == VTK_RECTILINEAR_GRID)
    {
        RectilinearGridExecute();
    }
    else if (do_type == VTK_STRUCTURED_GRID)
    {
        StructuredGridExecute();
    }
    else if (do_type == VTK_UNSTRUCTURED_GRID)
    {
        UnstructuredGridExecute();
    }
    else
    {
        GeneralExecute();
    }
110 111 112 113 114 115 116 117 118 119 120 121 122 123

    return 1;
}

// ****************************************************************************
//  Method: vtkSlicer::FillInputPortInformation
//
// ****************************************************************************

int
vtkSlicer::FillInputPortInformation(int, vtkInformation *info)
{
    info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
    return 1;
hrchilds's avatar
hrchilds committed
124 125
}

126
// ****************************************************************************
127
//  Class: SliceFunction
128
//
129 130
//  Purpose:
//    Slice functor that accesses points directly as array.
131
//
132
//  Notes:      
133
//
134 135 136 137
//  Programmer: Brad Whitlock
//  Creation:   Tue Mar 13 10:53:27 PDT 2012
//
//  Modifications:
138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177
//
// ****************************************************************************

template <typename T>
class SliceFunction
{
public:
    SliceFunction(const int *pt_dims, vtkPoints *pts, const double O[3], const double N[3])
    {
        pts_ptr = (const T *)pts->GetVoidPointer(0);
        ptstrideY = (vtkIdType)pt_dims[0];
        ptstrideZ = (vtkIdType)(pt_dims[0] * pt_dims[1]);
        Origin[0] = O[0]; Origin[1] = O[1]; Origin[2] = O[2];
        Normal[0] = N[0]; Normal[1] = N[1]; Normal[2] = N[2];
        D = Origin[0]*Normal[0] + Origin[1]*Normal[1] + Origin[2]*Normal[2];
    }

    inline T operator()(vtkIdType cellI,   vtkIdType cellJ,   vtkIdType cellK,
                        vtkIdType iOffset, vtkIdType jOffset, vtkIdType kOffset) const
    {
        vtkIdType ptId = (cellI + iOffset) + (cellJ + jOffset) * ptstrideY +
                         (cellK + kOffset) * ptstrideZ;
        const T *pt = pts_ptr + 3 * ptId;
        return pt[0]*Normal[0] + pt[1]*Normal[1] + pt[2]*Normal[2] - D;
    }

    inline T operator()(vtkIdType ptId) const
    {
        const T *pt = pts_ptr + 3 * ptId;
        return pt[0]*Normal[0] + pt[1]*Normal[1] + pt[2]*Normal[2] - D;
    }
private:
    vtkIdType ptstrideY, ptstrideZ; 
    const T *pts_ptr;
    double   Normal[3];
    double   Origin[3];
    double   D;
};

// ****************************************************************************
178
//  Class: GeneralSliceFunction
179
//
180 181
//  Purpose:
//    Slice functor that uses GetPoint to access points.
182
//
183
//  Notes:      
184
//
185 186 187 188
//  Programmer: Brad Whitlock
//  Creation:   Tue Mar 13 10:53:27 PDT 2012
//
//  Modifications:
189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209
//
// ****************************************************************************

class GeneralSliceFunction
{
public:
    GeneralSliceFunction(const int *pt_dims, vtkPoints *p, const double O[3], const double N[3])
    {
        pts = p;
        ptstrideY = (vtkIdType)pt_dims[0];
        ptstrideZ = (vtkIdType)(pt_dims[0] * pt_dims[1]);
        Origin[0] = O[0]; Origin[1] = O[1]; Origin[2] = O[2];
        Normal[0] = N[0]; Normal[1] = N[1]; Normal[2] = N[2];
        D = Origin[0]*Normal[0] + Origin[1]*Normal[1] + Origin[2]*Normal[2];
    }

    inline double operator()(vtkIdType cellI,   vtkIdType cellJ,   vtkIdType cellK,
                             vtkIdType iOffset, vtkIdType jOffset, vtkIdType kOffset) const
    {
        vtkIdType ptId = (cellI + iOffset) + (cellJ + jOffset) * ptstrideY +
                         (cellK + kOffset) * ptstrideZ;
210
        const double *pt = pts->GetPoint(ptId);
211 212 213 214 215
        return pt[0]*Normal[0] + pt[1]*Normal[1] + pt[2]*Normal[2] - D;
    }

    inline double operator()(vtkIdType ptId) const
    {
216
        const double *pt = pts->GetPoint(ptId);
217 218 219 220 221 222 223 224 225 226
        return pt[0]*Normal[0] + pt[1]*Normal[1] + pt[2]*Normal[2] - D;
    }
private:
    vtkIdType  ptstrideY, ptstrideZ; 
    vtkPoints *pts;
    double     Normal[3];
    double     Origin[3];
    double     D;
};

227 228 229 230 231 232
// ****************************************************************************
//  Method: vtkSlicer::StructuredGridExecute
//
//  Modifications:
//    Brad Whitlock, Thu Aug 12 14:51:27 PST 2004
//    Added float casts to the pow() arguments so it builds on MSVC7.Net.
hrchilds's avatar
hrchilds committed
233
//
234 235
//    Hank Childs, Sat Jan 27 12:45:03 PST 2007
//    Add check for 1xJxK and Ix1xK meshes (instead of crashing).
hrchilds's avatar
hrchilds committed
236
//
237 238 239
//    Brad Whitlock, Tue Mar 13 10:52:52 PDT 2012
//    I moved the implementation into vtkStructuredCreateTriangles and I added
//    different implementations.
240
//
241 242 243
//    Eric Brugger, Thu Jan 10 10:24:20 PST 2013
//    Modified to inherit from vtkPolyDataAlgorithm.
//
244 245 246
//    Brad Whitlock, Thu Jul 23 16:01:46 PDT 2015
//    Support for non-standard memory layout.
//
247
// ****************************************************************************
hrchilds's avatar
hrchilds committed
248 249 250 251

void
vtkSlicer::StructuredGridExecute(void)
{
252
    vtkStructuredGrid *sg = (vtkStructuredGrid *)input;
hrchilds's avatar
hrchilds committed
253 254
    int pt_dims[3];
    sg->GetDimensions(pt_dims);
hrchilds's avatar
hrchilds committed
255
    if (pt_dims[0] <= 1 || pt_dims[1] <= 1 || pt_dims[2] <= 1)
hrchilds's avatar
hrchilds committed
256 257 258 259 260
    {
        GeneralExecute();
        return;
    }

261
    vtkIdType          nCells = sg->GetNumberOfCells();
hrchilds's avatar
hrchilds committed
262 263 264 265
    vtkPoints         *inPts  = sg->GetPoints();
    vtkCellData       *inCD   = sg->GetCellData();
    vtkPointData      *inPD   = sg->GetPointData();

266
    vtkIdType ptSizeGuess = (this->CellList == NULL
hrchilds's avatar
hrchilds committed
267 268
                         ? (int) pow(float(nCells), 0.6667f) * 5 + 100
                         : CellListSize*5 + 100);
hrchilds's avatar
hrchilds committed
269 270 271

    vtkSurfaceFromVolume sfv(ptSizeGuess);

272 273 274 275 276 277 278 279 280 281
    int accessMethod = 0;
    if(inPts->GetData()->HasStandardMemoryLayout())
    {
        if(inPts->GetDataType() == VTK_FLOAT)
            accessMethod = 1;
        else if(inPts->GetDataType() == VTK_DOUBLE)
            accessMethod = 2;
    }

    if(accessMethod == 1)
hrchilds's avatar
hrchilds committed
282
    {
283 284 285 286 287
        vtkStructuredCreateTriangles<float, SliceFunction<float> >(
            sfv, this->CellList, this->CellListSize, nCells,
            pt_dims, SliceFunction<float>(pt_dims, inPts, this->Origin, this->Normal)
        );
    }
288
    else if(accessMethod == 2)
289 290 291 292 293 294 295 296 297 298 299 300
    {
        vtkStructuredCreateTriangles<double, SliceFunction<double> >(
            sfv, this->CellList, this->CellListSize, nCells,
            pt_dims, SliceFunction<double>(pt_dims, inPts, this->Origin, this->Normal)
        );
    }
    else
    {
        vtkStructuredCreateTriangles<double, GeneralSliceFunction >(
            sfv, this->CellList, this->CellListSize, nCells,
            pt_dims, GeneralSliceFunction(pt_dims, inPts, this->Origin, this->Normal)
        );
hrchilds's avatar
hrchilds committed
301 302
    }

303
    sfv.ConstructPolyData(inPD, inCD, output, inPts);
hrchilds's avatar
hrchilds committed
304 305
}

306
// ****************************************************************************
307
//  Class: RectSliceFunction
308
//
309 310
//  Purpose:
//    Rectilinear slice function that accesses the coordinate arrays directly.
311
//
312
//  Notes:      
313
//
314 315 316 317
//  Programmer: Brad Whitlock
//  Creation:   Tue Mar 13 10:50:45 PDT 2012
//
//  Modifications:
318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361
//
// ****************************************************************************

template <typename T>
class RectSliceFunction
{
public:
    RectSliceFunction(vtkDataArray *Xc, vtkDataArray *Yc, 
        vtkDataArray *Zc, const double O[3], const double N[3])
    {
        X = NULL;
        if(Xc != NULL)
            X = (const T *)Xc->GetVoidPointer(0);
        Y = NULL;
        if(Yc != NULL)
            Y = (const T *)Yc->GetVoidPointer(0);
        Z = NULL;
        if(Zc != NULL)
            Z = (const T *)Zc->GetVoidPointer(0);

        Origin[0] = O[0]; Origin[1] = O[1]; Origin[2] = O[2];
        Normal[0] = N[0]; Normal[1] = N[1]; Normal[2] = N[2];
        D = Origin[0]*Normal[0] + Origin[1]*Normal[1] + Origin[2]*Normal[2];
    }

    inline T operator()(vtkIdType cellI,   vtkIdType cellJ,   vtkIdType cellK,
                        vtkIdType iOffset, vtkIdType jOffset, vtkIdType kOffset) const
    {
        T pt[3];
        pt[0] = X[cellI + iOffset];
        pt[1] = Y[cellJ + jOffset];
        pt[2] = Z[cellK + kOffset];
        return pt[0]*Normal[0] + pt[1]*Normal[1] + pt[2]*Normal[2] - D;
    }
private:
    const T *X;
    const T *Y;
    const T *Z;
    double   Normal[3];
    double   Origin[3];
    double   D;
};

// ****************************************************************************
362
//  Class: GeneralRectSliceFunction
363
//
364 365 366
//  Purpose:
//    Rectilinear slice function that uses GetTuple to access the coordinate
//     arrays.
367
//
368
//  Notes:      
369
//
370 371 372 373
//  Programmer: Brad Whitlock
//  Creation:   Tue Mar 13 10:50:19 PDT 2012
//
//  Modifications:
374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408
//
// ****************************************************************************

class GeneralRectSliceFunction
{
public:
    GeneralRectSliceFunction(vtkDataArray *Xc, vtkDataArray *Yc, 
        vtkDataArray *Zc, const double O[3], const double N[3])
    {
        X = Xc;
        Y = Yc;
        Z = Zc;
        Origin[0] = O[0]; Origin[1] = O[1]; Origin[2] = O[2];
        Normal[0] = N[0]; Normal[1] = N[1]; Normal[2] = N[2];
        D = Origin[0]*Normal[0] + Origin[1]*Normal[1] + Origin[2]*Normal[2];
    }

    inline double operator()(vtkIdType cellI,   vtkIdType cellJ,   vtkIdType cellK,
                             vtkIdType iOffset, vtkIdType jOffset, vtkIdType kOffset) const
    {
        double pt[3];
        pt[0] = X->GetTuple1(cellI + iOffset);
        pt[1] = Y->GetTuple1(cellJ + jOffset);
        pt[2] = Z->GetTuple1(cellK + kOffset);
        return pt[0]*Normal[0] + pt[1]*Normal[1] + pt[2]*Normal[2] - D;
    }
private:
    vtkDataArray *X;
    vtkDataArray *Y;
    vtkDataArray *Z;
    double        Normal[3];
    double        Origin[3];
    double        D;
};

409 410 411 412 413 414 415 416 417
// ****************************************************************************
//  Method: vtkSlicer::RectilinearGridExecute
//
//  Modifications:
//    Brad Whitlock, Thu Aug 12 14:51:27 PST 2004
//    Added float casts to the pow() arguments so it builds on MSVC7.Net.
//
//    Hank Childs, Sat Jan 27 12:45:03 PST 2007
//    Add check for 1xJxK and Ix1xK meshes (instead of crashing).
hrchilds's avatar
hrchilds committed
418
//
419 420 421 422
//    Brad Whitlock, Tue Mar 13 10:48:45 PDT 2012
//    Moved the implementation into vtkStructuredCreateTriangles and added
//    different implementations for float/double and a general method that
//    uses GetTuple to get the values.
hrchilds's avatar
hrchilds committed
423
//
424 425
//    Eric Brugger, Thu Jan 10 10:24:20 PST 2013
//    Modified to inherit from vtkPolyDataAlgorithm.
426
//
427 428 429
//    Brad Whitlock, Thu Jul 23 16:01:46 PDT 2015
//    Support for non-standard memory layout.
//
430
// ****************************************************************************
hrchilds's avatar
hrchilds committed
431

432 433
void
vtkSlicer::RectilinearGridExecute(void)
hrchilds's avatar
hrchilds committed
434
{
435
    vtkRectilinearGrid *rg = (vtkRectilinearGrid *)input;
hrchilds's avatar
hrchilds committed
436 437
    int pt_dims[3];
    rg->GetDimensions(pt_dims);
hrchilds's avatar
hrchilds committed
438
    if (pt_dims[0] <= 1 || pt_dims[1] <= 1 || pt_dims[2] <= 1)
hrchilds's avatar
hrchilds committed
439 440 441 442 443
    {
        GeneralExecute();
        return;
    }

444
    vtkIdType     nCells = rg->GetNumberOfCells();
hrchilds's avatar
hrchilds committed
445 446 447
    vtkCellData  *inCD   = rg->GetCellData();
    vtkPointData *inPD   = rg->GetPointData();

448
    vtkIdType ptSizeGuess = (this->CellList == NULL
hrchilds's avatar
hrchilds committed
449 450
                         ? (int) pow(float(nCells), 0.6667f) * 5 + 100
                         : CellListSize*5 + 100);
hrchilds's avatar
hrchilds committed
451 452 453

    vtkSurfaceFromVolume sfv(ptSizeGuess);

454 455 456
    int tx = rg->GetXCoordinates()->GetDataType();
    int ty = rg->GetYCoordinates()->GetDataType();
    int tz = rg->GetZCoordinates()->GetDataType();
457 458 459 460 461 462
    bool smlx = rg->GetXCoordinates()->HasStandardMemoryLayout();
    bool smly = rg->GetYCoordinates()->HasStandardMemoryLayout();
    bool smlz = rg->GetZCoordinates()->HasStandardMemoryLayout();
    bool sameTypes = tx == ty && ty == tz;
    bool sameML = smlx == smly && smly == smlz;
    bool same = sameTypes && sameML;
hrchilds's avatar
hrchilds committed
463

464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489
    if(same && tx == VTK_FLOAT)
    {
        vtkStructuredCreateTriangles<float, RectSliceFunction<float> >(
            sfv, this->CellList, this->CellListSize, nCells, pt_dims,
            RectSliceFunction<float>(
                rg->GetXCoordinates(), rg->GetYCoordinates(), rg->GetZCoordinates(),
                this->Origin, this->Normal)
        );
    }
    else if(same && tx == VTK_DOUBLE)
    {
        vtkStructuredCreateTriangles<double, RectSliceFunction<double> >(
            sfv, this->CellList, this->CellListSize, nCells, pt_dims,
            RectSliceFunction<double>(
                rg->GetXCoordinates(), rg->GetYCoordinates(), rg->GetZCoordinates(),
                this->Origin, this->Normal)
        );
    }
    else
    {
        vtkStructuredCreateTriangles<double, GeneralRectSliceFunction >(
            sfv, this->CellList, this->CellListSize, nCells, pt_dims,
            GeneralRectSliceFunction(
                rg->GetXCoordinates(), rg->GetYCoordinates(), rg->GetZCoordinates(),
                this->Origin, this->Normal)
        );
hrchilds's avatar
hrchilds committed
490 491
    }

492 493
    sfv.ConstructPolyData(inPD, inCD, output, pt_dims, 
        rg->GetXCoordinates(), rg->GetYCoordinates(), rg->GetZCoordinates());
hrchilds's avatar
hrchilds committed
494 495
}

hrchilds's avatar
hrchilds committed
496
// ****************************************************************************
497
//  Method: vtkSlicer::UnstructuredGridExecute
hrchilds's avatar
hrchilds committed
498
//
499
//  Modifications:
hrchilds's avatar
hrchilds committed
500 501 502
//    Hank Childs, Tue Mar 30 07:07:42 PST 2004
//    Add support for slicing vertices.
//
hrchilds's avatar
hrchilds committed
503 504 505
//    Brad Whitlock, Thu Aug 12 14:51:27 PST 2004
//    Added float casts to the pow() arguments so it builds on MSVC7.Net.
//
506 507 508
//    Brad Whitlock, Tue Mar 13 16:55:36 PDT 2012
//    I moved code to generate triangles into vtkUnstructuredCreateTriangles.
//
509 510 511
//    Brad Whitlock, Wed Apr 11 11:37:18 PDT 2012
//    When we can't slice a cell, insert faces too for polyhedral cells.
//
512 513 514
//    Eric Brugger, Thu Jan 10 10:24:20 PST 2013
//    Modified to inherit from vtkPolyDataAlgorithm.
//
515 516 517
//    Kathleen Biagas, Fri Jan 25 16:04:46 PST 2013
//    Call Update on the filter, not the data object.
//
hrchilds's avatar
hrchilds committed
518 519
// ****************************************************************************

520 521
void
vtkSlicer::UnstructuredGridExecute(void)
hrchilds's avatar
hrchilds committed
522 523 524 525 526 527 528 529 530 531 532
{
    // The routine here is a bit trickier than for the Rectilinear or
    // Structured grids.  We want to slice an unstructured grid -- but that
    // could mean any cell type.  We only have triangulation tables for
    // the finite element zoo.  So the gameplan is to slice any of the
    // elements of the finite element zoo.  If there are more elements left
    // over, slice them using the conventional VTK filters.  Finally,
    // append together the slices from the zoo with the slices from the
    // non-zoo elements.  If all the elements are from the zoo, then just
    // slice them with no appending.

533
    vtkUnstructuredGrid *ug = (vtkUnstructuredGrid *)input;
hrchilds's avatar
hrchilds committed
534

535
    vtkIdType          nCells = ug->GetNumberOfCells();
hrchilds's avatar
hrchilds committed
536 537 538 539
    vtkPoints         *inPts  = ug->GetPoints();
    vtkCellData       *inCD   = ug->GetCellData();
    vtkPointData      *inPD   = ug->GetPointData();

540
    vtkIdType ptSizeGuess = (this->CellList == NULL
hrchilds's avatar
hrchilds committed
541 542
                         ? (int) pow(float(nCells), 0.6667f) * 5 + 100
                         : CellListSize*5 + 100);
hrchilds's avatar
hrchilds committed
543 544 545 546

    vtkSurfaceFromVolume sfv(ptSizeGuess);

    vtkUnstructuredGrid *stuff_I_cant_slice = vtkUnstructuredGrid::New();
hrchilds's avatar
hrchilds committed
547
    vtkPolyData *vertices_on_slice = vtkPolyData::New();
hrchilds's avatar
hrchilds committed
548

hrchilds's avatar
hrchilds committed
549
    double D = Origin[0]*Normal[0] + Origin[1]*Normal[1] + Origin[2]*Normal[2];
hrchilds's avatar
hrchilds committed
550

551 552 553 554
    vtkIdType nToProcess = (CellList != NULL ? CellListSize : nCells);
    vtkIdType numIcantSlice = 0;
    vtkIdType numVertices = 0;
    for (vtkIdType i = 0 ; i < nToProcess ; i++)
hrchilds's avatar
hrchilds committed
555
    {
556
        vtkIdType  cellId = (CellList != NULL ? CellList[i] : i);
hrchilds's avatar
hrchilds committed
557
        int        cellType = ug->GetCellType(cellId);
558
        vtkIdType  npts;
hrchilds's avatar
hrchilds committed
559 560
        vtkIdType *pts;
        ug->GetCellPoints(cellId, npts, pts);
561 562
        const int *triangulation_table = NULL;
        const int *vertices_from_edges = NULL;
hrchilds's avatar
hrchilds committed
563 564
        int tt_step = 0;
        bool canSlice = false;
hrchilds's avatar
hrchilds committed
565
        bool isVertex = false;
hrchilds's avatar
hrchilds committed
566 567 568
        switch (cellType)
        {
          case VTK_TETRA:
569 570
            triangulation_table = (const int *) tetTriangulationTable;
            vertices_from_edges = (const int *) tetVerticesFromEdges;
hrchilds's avatar
hrchilds committed
571 572 573 574 575
            tt_step = 7;
            canSlice = true;
            break;
 
          case VTK_PYRAMID:
576 577
            triangulation_table = (const int *) pyramidTriangulationTable;
            vertices_from_edges = (const int *) pyramidVerticesFromEdges;
hrchilds's avatar
hrchilds committed
578 579 580 581 582
            tt_step = 13;
            canSlice = true;
            break;
 
          case VTK_WEDGE:
583 584
            triangulation_table = (const int *) wedgeTriangulationTable;
            vertices_from_edges = (const int *) wedgeVerticesFromEdges;
hrchilds's avatar
hrchilds committed
585 586 587 588 589
            tt_step = 13;
            canSlice = true;
            break;
 
          case VTK_HEXAHEDRON:
590 591
            triangulation_table = (const int *) hexTriangulationTable;
            vertices_from_edges = (const int *) hexVerticesFromEdges;
hrchilds's avatar
hrchilds committed
592 593 594
            tt_step = 16;
            canSlice = true;
            break;
hrchilds's avatar
hrchilds committed
595 596 597 598
 
          case VTK_VERTEX:
            isVertex = true;
            break;
hrchilds's avatar
hrchilds committed
599 600 601 602 603 604 605 606

          default:
            canSlice = false;
            break;
        }
 
        if (canSlice)
        {
607 608
            int tmp[3] = {0,0,0};
            if(inPts->GetDataType() == VTK_FLOAT)
hrchilds's avatar
hrchilds committed
609
            {
610 611 612 613 614
                vtkUnstructuredCreateTriangles<float, SliceFunction<float> >(
                    sfv, cellId, pts, npts,
                    triangulation_table, vertices_from_edges, tt_step,
                    SliceFunction<float>(tmp, inPts, this->Origin, this->Normal)
                );
hrchilds's avatar
hrchilds committed
615
            }
616
            else if(inPts->GetDataType() == VTK_DOUBLE)
hrchilds's avatar
hrchilds committed
617
            {
618 619 620 621 622 623 624 625 626 627 628 629 630
                vtkUnstructuredCreateTriangles<double, SliceFunction<double> >(
                    sfv, cellId, pts, npts,
                    triangulation_table, vertices_from_edges, tt_step,
                    SliceFunction<double>(tmp, inPts, this->Origin, this->Normal)
                );
            }
            else
            {
                vtkUnstructuredCreateTriangles<double, GeneralSliceFunction >(
                    sfv, cellId, pts, npts,
                    triangulation_table, vertices_from_edges, tt_step,
                    GeneralSliceFunction(tmp, inPts, this->Origin, this->Normal)
                );
hrchilds's avatar
hrchilds committed
631 632
            }
        }
hrchilds's avatar
hrchilds committed
633 634 635 636 637
        else if (isVertex)
        {
            //
            // Determine if the vertex is even on the plane.
            //
638
            const double *pt = inPts->GetPoint(pts[0]);
hrchilds's avatar
hrchilds committed
639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657
            double dist_from_plane = Normal[0]*pt[0] + Normal[1]*pt[1]
                                   + Normal[2]*pt[2] - D;
            if (fabs(dist_from_plane) < 1e-12)
            {
                if (numVertices == 0)
                {
                    vertices_on_slice->SetPoints(ug->GetPoints());
                    vertices_on_slice->GetPointData()->ShallowCopy(
                                                           ug->GetPointData());
                    vertices_on_slice->Allocate(nCells);
                    vertices_on_slice->GetCellData()->
                                       CopyAllocate(ug->GetCellData(), nCells);
                }
                vertices_on_slice->InsertNextCell(VTK_VERTEX, 1, pts);
                vertices_on_slice->GetCellData()->
                              CopyData(ug->GetCellData(), cellId, numVertices);
                numVertices++;
            }
        }
hrchilds's avatar
hrchilds committed
658 659 660
        else
        {
            if (numIcantSlice == 0)
hrchilds's avatar
hrchilds committed
661 662 663 664 665
            {
                stuff_I_cant_slice->SetPoints(ug->GetPoints());
                stuff_I_cant_slice->GetPointData()->ShallowCopy(
                                                           ug->GetPointData());
                stuff_I_cant_slice->Allocate(nCells);
hrchilds's avatar
hrchilds committed
666 667
                stuff_I_cant_slice->GetCellData()->
                                       CopyAllocate(ug->GetCellData(), nCells);
hrchilds's avatar
hrchilds committed
668
            }
hrchilds's avatar
hrchilds committed
669

670 671 672 673 674 675 676 677 678
            if(cellType == VTK_POLYHEDRON)
            {
                vtkIdType nFaces, *facePtIds;
                ug->GetFaceStream(cellId, nFaces, facePtIds);
                stuff_I_cant_slice->InsertNextCell(cellType, npts, pts, 
                    nFaces, facePtIds);
            }
            else
                stuff_I_cant_slice->InsertNextCell(cellType, npts, pts);
hrchilds's avatar
hrchilds committed
679 680 681 682 683 684
            stuff_I_cant_slice->GetCellData()->
                            CopyData(ug->GetCellData(), cellId, numIcantSlice);
            numIcantSlice++;
        }
    }

hrchilds's avatar
hrchilds committed
685
    if ((numIcantSlice > 0) || (numVertices > 0))
hrchilds's avatar
hrchilds committed
686
    {
hrchilds's avatar
hrchilds committed
687 688 689 690 691
        vtkAppendPolyData *appender = vtkAppendPolyData::New();

        if (numIcantSlice > 0)
        {
            vtkPolyData *not_from_zoo  = vtkPolyData::New();
692
            SliceDataset(stuff_I_cant_slice, not_from_zoo, true);
693
            appender->AddInputData(not_from_zoo);
hrchilds's avatar
hrchilds committed
694 695 696 697 698
            not_from_zoo->Delete();
        }

        if (numVertices > 0)
        {
699
            appender->AddInputData(vertices_on_slice);
hrchilds's avatar
hrchilds committed
700 701
        }

hrchilds's avatar
hrchilds committed
702
        vtkPolyData *just_from_zoo = vtkPolyData::New();
703
        sfv.ConstructPolyData(inPD, inCD, just_from_zoo, inPts);
704
        appender->AddInputData(just_from_zoo);
hrchilds's avatar
hrchilds committed
705 706
        just_from_zoo->Delete();

707
        appender->Update();
hrchilds's avatar
hrchilds committed
708

hrchilds's avatar
hrchilds committed
709
        output->ShallowCopy(appender->GetOutput());
hrchilds's avatar
hrchilds committed
710 711 712 713
        appender->Delete();
    }
    else
    {
714
        sfv.ConstructPolyData(inPD, inCD, output, inPts);
hrchilds's avatar
hrchilds committed
715 716 717
    }

    stuff_I_cant_slice->Delete();
hrchilds's avatar
hrchilds committed
718
    vertices_on_slice->Delete();
hrchilds's avatar
hrchilds committed
719 720
}

721 722 723 724 725 726
// ****************************************************************************
//  Modifications:
//    Eric Brugger, Thu Jan 10 10:24:20 PST 2013
//    Modified to inherit from vtkPolyDataAlgorithm.
//
// ****************************************************************************
hrchilds's avatar
hrchilds committed
727

728 729
void
vtkSlicer::GeneralExecute(void)
hrchilds's avatar
hrchilds committed
730
{
731
    SliceDataset(input, output, false);
hrchilds's avatar
hrchilds committed
732 733
}

hrchilds's avatar
hrchilds committed
734 735 736 737 738
// ****************************************************************************
//  Modifications:
//    Kathleen Bonnell, Wed Apr 27 18:47:18 PDT 2005
//    Use vtkVisItCutter, which has modifications to correctly handle CellData. 
//
739 740 741 742
//    Brad Whitlock, Wed Apr 11 11:34:16 PDT 2012
//    Use vtkCutter when we're taking this path so we can slice general VTK
//    datasets containing cells that vtkVisItCutter can't slice.
//
hrchilds's avatar
hrchilds committed
743
// ****************************************************************************
744

745 746
void
vtkSlicer::SliceDataset(vtkDataSet *in_ds, vtkPolyData *out_pd,
747
    bool useVTKFilter)
hrchilds's avatar
hrchilds committed
748 749 750 751 752
{
    vtkPlane  *plane  = vtkPlane::New();
    plane->SetOrigin(Origin[0], Origin[1], Origin[2]);
    plane->SetNormal(Normal[0], Normal[1], Normal[2]);

753 754 755 756
    if(useVTKFilter)
    {
        vtkCutter *cutter = vtkCutter::New();
        cutter->SetCutFunction(plane);
757
        cutter->SetInputData(in_ds);
758 759 760 761 762 763 764 765 766
        cutter->Update();

        out_pd->ShallowCopy(cutter->GetOutput());
        cutter->Delete();
    }
    else
    {
        vtkVisItCutter *cutter = vtkVisItCutter::New();
        cutter->SetCutFunction(plane);
767
        cutter->SetInputData(in_ds);
768 769 770 771 772
        cutter->Update();

        out_pd->ShallowCopy(cutter->GetOutput());
        cutter->Delete();
    }
hrchilds's avatar
hrchilds committed
773

hrchilds's avatar
hrchilds committed
774 775 776
    plane->Delete();
}

777 778 779 780 781 782 783
// ****************************************************************************
//  Method: vtkSlicer::PrintSelf
//
// ****************************************************************************

void
vtkSlicer::PrintSelf(ostream& os, vtkIndent indent)
hrchilds's avatar
hrchilds committed
784 785 786 787 788 789 790 791
{
  this->Superclass::PrintSelf(os,indent);

  os << indent << "Normal: " << Normal[0] << ", " << Normal[1] << ", "
     << Normal[2] << "\n";
  os << indent << "Origin: " << Origin[0] << ", " << Origin[1] << ", "
     << Origin[2] << "\n";
}