vtkVisItClipper.C 32 KB
Newer Older
hrchilds's avatar
hrchilds committed
1 2
/*****************************************************************************
*
3
* Copyright (c) 2000 - 2009, Lawrence Livermore National Security, LLC
hrchilds's avatar
hrchilds committed
4
* Produced at the Lawrence Livermore National Laboratory
5
* LLNL-CODE-400124
hrchilds's avatar
hrchilds committed
6 7
* All rights reserved.
*
8
* This file is  part of VisIt. For  details, see https://visit.llnl.gov/.  The
hrchilds's avatar
hrchilds committed
9 10 11 12 13 14 15 16 17 18
* full copyright notice is contained in the file COPYRIGHT located at the root
* of the VisIt distribution or at http://www.llnl.gov/visit/copyright.html.
*
* Redistribution  and  use  in  source  and  binary  forms,  with  or  without
* modification, are permitted provided that the following conditions are met:
*
*  - Redistributions of  source code must  retain the above  copyright notice,
*    this list of conditions and the disclaimer below.
*  - Redistributions in binary form must reproduce the above copyright notice,
*    this  list of  conditions  and  the  disclaimer (as noted below)  in  the
19 20 21
*    documentation and/or other materials provided with the distribution.
*  - Neither the name of  the LLNS/LLNL nor the names of  its contributors may
*    be used to endorse or promote products derived from this software without
hrchilds's avatar
hrchilds committed
22 23 24 25 26
*    specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT  HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR  PURPOSE
27 28 29
* ARE  DISCLAIMED. IN  NO EVENT  SHALL LAWRENCE  LIVERMORE NATIONAL  SECURITY,
* LLC, THE  U.S.  DEPARTMENT OF  ENERGY  OR  CONTRIBUTORS BE  LIABLE  FOR  ANY
* DIRECT,  INDIRECT,   INCIDENTAL,   SPECIAL,   EXEMPLARY,  OR   CONSEQUENTIAL
hrchilds's avatar
hrchilds committed
30 31 32 33 34 35 36 37 38
* DAMAGES (INCLUDING, BUT NOT  LIMITED TO, PROCUREMENT OF  SUBSTITUTE GOODS OR
* SERVICES; LOSS OF  USE, DATA, OR PROFITS; OR  BUSINESS INTERRUPTION) HOWEVER
* CAUSED  AND  ON  ANY  THEORY  OF  LIABILITY,  WHETHER  IN  CONTRACT,  STRICT
* LIABILITY, OR TORT  (INCLUDING NEGLIGENCE OR OTHERWISE)  ARISING IN ANY  WAY
* OUT OF THE  USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
* DAMAGE.
*
*****************************************************************************/

hrchilds's avatar
hrchilds committed
39
#include "vtkVisItClipper.h"
hrchilds's avatar
hrchilds committed
40
#include <vtkAppendFilter.h>
hrchilds's avatar
hrchilds committed
41 42 43
#include <vtkCellData.h>
#include <vtkClipDataSet.h>
#include <vtkFloatArray.h>
hrchilds's avatar
hrchilds committed
44
#include <vtkImplicitFunction.h>
hrchilds's avatar
hrchilds committed
45
#include <vtkObjectFactory.h>
hrchilds's avatar
hrchilds committed
46
#include <vtkPlane.h>
hrchilds's avatar
hrchilds committed
47
#include <vtkPointData.h>
hrchilds's avatar
hrchilds committed
48
#include <vtkPolyData.h>
hrchilds's avatar
hrchilds committed
49
#include <vtkQuadric.h>
hrchilds's avatar
hrchilds committed
50 51 52
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>
#include <vtkUnstructuredGrid.h>
hrchilds's avatar
hrchilds committed
53 54
#include <vtkVolumeFromVolume.h>

hrchilds's avatar
hrchilds committed
55 56 57 58 59 60 61
#include <ImproperUseException.h>

#include <DebugStream.h>

#include <math.h>
#include <vector>

hrchilds's avatar
hrchilds committed
62
#include <ClipCases.h>
hrchilds's avatar
hrchilds committed
63 64
#include <vtkTriangulationTables.h>

hrchilds's avatar
hrchilds committed
65 66
vtkCxxRevisionMacro(vtkVisItClipper, "$Revision: 1.00 $");
vtkStandardNewMacro(vtkVisItClipper);
hrchilds's avatar
hrchilds committed
67

hrchilds's avatar
hrchilds committed
68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128
//
// Function: AdjustPercentToZeroCrossing
//
// Purpose: Given coordinate array, point ids and linear estimate of
// a cut, use quadric to compute actual zero crossing and adjust the
// percent value to hit the zero crossing
//
//  Programmer: Mark C. Miller
//  Creation:   December 3, 2006 
//
static void
AdjustPercentToZeroCrossing(const float *const pts, int ptId1, int ptId2,
    vtkImplicitFunction *func, float *percent)
{
    if (func == 0)
        return;

    // we only handle general quadrics at the moment
    if (strcmp(func->GetClassName(), "vtkQuadric") != 0)
        return;

    //
    // quadric equation coefficient array indexing...
    // x^2   y^2   z^2    xy    xz    yz    x    y    z    1
    //  0     1     2     3     4     5     6    7    8    9
    //
    vtkQuadric *quadric = vtkQuadric::SafeDownCast(func);
    const double *a = quadric->GetCoefficients();

    // quick check for planar functions. They're linear and so
    // 'percent' is already correct
    if (a[0] == 0.0 && a[1] == 0.0 && a[2] == 0.0 &&
        a[3] == 0.0 && a[4] == 0.0 && a[5] == 0.0)
        return;

    //
    // We'll define a "ray" between points p0 and p1 such that a
    // point along it is defined by p(t) = p0 + t * (p1 - p0).
    // When t==0, p(t)==p0 and when t==1, p(t)==p1. So, along
    // the edge between the points p0 and p1, 0<=t<=1
    //
    const float *const p0 = pts + 3*ptId1;
    const float *const p1 = pts + 3*ptId2;

    // origin of "ray" to intersect against the quadric surface
    double x0 = p0[0];
    double y0 = p0[1];
    double z0 = p0[2];

    // direction (non-normalized) of ray to intersect quadric surface
    double xd = p1[0] - x0;
    double yd = p1[1] - y0;
    double zd = p1[2] - z0;

    //
    // compute quadratic equation coefficients for ray/quadric intersection
    // At^2 + Bt + C = 0
    //
    // These equations were obtained from various web resources. However,
    // I am suspect of the equation for the B coefficient as cited on the
    // web. Several sources cite the equation with the commented line. However,
ahern's avatar
ahern committed
129
    // there is an asymmetry in it where the coefficient of the a[5] term does
hrchilds's avatar
hrchilds committed
130 131 132 133 134 135 136 137 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
    // not include a xd*z0 contribution analagous to the a[3] and a[4] terms.
    // Empirical results from its use have shown that indeed it is in error.
    // The commented line and this comment is left here in case anyone
    // bothers to check this math against available sources.
    //
    double A = a[0]*xd*xd + a[1]*yd*yd + a[2]*zd*zd +
               a[3]*xd*yd + a[4]*yd*zd + a[5]*xd*zd;
    double B = 2*a[0]*x0*xd + 2*a[1]*y0*yd + 2*a[2]*z0*zd +
               //a[3]*(x0*yd+y0*xd) + a[4]*(y0*zd+yd*z0) + a[5]*x0*zd +
               a[3]*(x0*yd+xd*y0) + a[4]*(y0*zd+yd*z0) + a[5]*(x0*zd+xd*z0) +
               a[6]*xd + a[7]*yd +a[8]*zd;
    double C = a[0]*x0*x0 + a[1]*y0*y0 + a[2]*z0*z0 +
               a[3]*x0*y0 + a[4]*y0*z0 + a[5]*x0*z0 +
               a[6]*x0 + a[7]*y0 + a[8]*z0 + a[9];

    //
    // compute the root(s) of the quadratic equation
    //
    double t = 0.0;
    if (A == 0)
    {
        //
        // We get here if the quadric is really just linear
        //
        if (B == 0)
            t = 0.0;
        else
            t = -C / B;
    }
    else
    {
        //
        // We get here only when the quadric is indeed non-linear
        //
        double disc = B*B - 4*A*C;
        if (disc >= 0.0)
        {
            t = (-B - sqrt(disc)) / (2*A);
            if (t < 0)
                t = (-B + sqrt(disc)) / (2*A);
        }
    }

    if (t > 0.0 && t <= 1.0)
        *percent = 1.0-t;
}

hrchilds's avatar
hrchilds committed
177 178 179 180 181 182 183 184 185 186
// ****************************************************************************
//  Constructor:  vtkVisItClipper::vtkVisItClipper
//
//  Programmer:  Jeremy Meredith
//  Creation:    August 11, 2003
//
//  Modifications:
//    Jeremy Meredith, Tue Aug 29 13:38:08 EDT 2006
//    Added support for leaving cells whole.
//
187 188 189
//    Hank Childs, Sat Sep 29 11:14:58 PDT 2007
//    Initialize new data members.
//
hrchilds's avatar
hrchilds committed
190
// ****************************************************************************
191

hrchilds's avatar
hrchilds committed
192
vtkVisItClipper::vtkVisItClipper()
hrchilds's avatar
hrchilds committed
193 194 195 196 197
{
    CellList = NULL;
    CellListSize = 0;
    insideOut = false;
    clipFunction = NULL;
hrchilds's avatar
hrchilds committed
198
    removeWholeCells = false;
hrchilds's avatar
hrchilds committed
199 200 201
    useZeroCrossings = false;
    computeInsideAndOut = false;
    otherOutput = NULL;
202 203
    scalarArrayAsVTK = NULL;
    iOwnData = false;
hrchilds's avatar
hrchilds committed
204 205
}

hrchilds's avatar
hrchilds committed
206 207 208 209 210 211 212 213
// ****************************************************************************
//  Destructor:  vtkVisItClipper::~vtkVisItClipper
//
//  Programmer:  Jeremy Meredith
//  Creation:    August 11, 2003
//
//  Modifications:
//
214 215 216
//    Hank Childs, Sat Sep 29 11:14:58 PDT 2007
//    Clean up new data members.
//
hrchilds's avatar
hrchilds committed
217
// ****************************************************************************
hrchilds's avatar
hrchilds committed
218
vtkVisItClipper::~vtkVisItClipper()
hrchilds's avatar
hrchilds committed
219
{
hrchilds's avatar
hrchilds committed
220 221
    if (otherOutput)
        otherOutput->Delete();
222 223 224 225
    if (iOwnData)
        delete [] scalarArray;
    if (scalarArrayAsVTK != NULL)
        scalarArrayAsVTK->Delete();
hrchilds's avatar
hrchilds committed
226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245
}

void
vtkVisItClipper::SetUseZeroCrossings(bool use)
{
    if (use && clipFunction && 
        (strcmp(clipFunction->GetClassName(), "vtkQuadric") != 0))
    {
        vtkErrorMacro("UseZeroCrossings set to true allowed only with "
                      "vtkQuadric implicit functions");
        return;
    }

    useZeroCrossings = use;
}

void
vtkVisItClipper::SetComputeInsideAndOut(bool compute)
{
    computeInsideAndOut = compute;
hrchilds's avatar
hrchilds committed
246 247 248
}

void
hrchilds's avatar
hrchilds committed
249
vtkVisItClipper::SetCellList(int *cl, int size)
hrchilds's avatar
hrchilds committed
250 251 252 253 254 255
{
    CellList = cl;
    CellListSize = size;
}

void
hrchilds's avatar
hrchilds committed
256
vtkVisItClipper::SetClipFunction(vtkImplicitFunction *func)
hrchilds's avatar
hrchilds committed
257
{
hrchilds's avatar
hrchilds committed
258 259 260 261 262 263 264
    if (useZeroCrossings && (strcmp(func->GetClassName(), "vtkQuadric") != 0))
    {
        vtkErrorMacro("Only vtkQuadric implicit functions "
                      "allowed with UseZeroCrossings set to true");
        return;
    }

hrchilds's avatar
hrchilds committed
265
    // Set the clip function
hrchilds's avatar
hrchilds committed
266
    clipFunction = func;
hrchilds's avatar
hrchilds committed
267 268 269 270 271

    // Clear the scalar array so we know to use the clip function
    scalarArray = NULL;
}

hrchilds's avatar
hrchilds committed
272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287
// ****************************************************************************
//  Method:  vtkVisItClipper::SetClipScalars
//
//  Purpose:
//    Set the scalar array used for clipping, and the cutoff.
//    To clip to a range, execute this filter once for the minimum
//    and once for the maximum.
//
//  Arguments:
//    array      the scalar array
//    cutoff     the cutoff
//
//  Programmer:  Jeremy Meredith
//  Creation:    January 30, 2004
//
//  Modifications:
288
//
hrchilds's avatar
hrchilds committed
289 290 291 292
//    Jeremy Meredith, Wed May  5 14:48:23 PDT 2004
//    Made it allow only a single cutoff, and use the "insideOut"
//    value to determine if this is a min or max value.
//
293 294 295 296
//    Hank Childs, Sat Sep 29 11:14:58 PDT 2007
//    Change the array argument to be a vtk data type.  Also added support
//    for data types besides "float".
//
hrchilds's avatar
hrchilds committed
297
// ****************************************************************************
298

hrchilds's avatar
hrchilds committed
299
void
300
vtkVisItClipper::SetClipScalars(vtkDataArray *array, float cutoff)
hrchilds's avatar
hrchilds committed
301
{
302 303 304 305 306 307 308 309 310 311 312
    if (iOwnData)
    {
        delete [] scalarArray;
        iOwnData = false;
    }
    if (scalarArrayAsVTK != NULL)
    {
        scalarArrayAsVTK->Delete();
        scalarArrayAsVTK = NULL;
    }

hrchilds's avatar
hrchilds committed
313 314 315 316
    // Clear the clip function so we know to use scalars
    clipFunction = NULL;

    // Set the scalar array
317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332
    scalarArrayAsVTK = array;
    scalarArrayAsVTK->Register(NULL);
    if (array->GetDataType() == VTK_FLOAT)
    {
        scalarArray = (float *) array->GetVoidPointer(0);
    }
    else
    {
        iOwnData = true;
        int nTuples = array->GetNumberOfTuples();
        scalarArray = new float[nTuples];
        for (int i = 0 ; i < nTuples ; i++)
        {
            scalarArray[i] = array->GetTuple1(i);
        }
    }
hrchilds's avatar
hrchilds committed
333

hrchilds's avatar
hrchilds committed
334 335
    // Set the cutoff
    scalarCutoff     = cutoff;
hrchilds's avatar
hrchilds committed
336 337 338
}

void
hrchilds's avatar
hrchilds committed
339
vtkVisItClipper::SetInsideOut(bool io)
hrchilds's avatar
hrchilds committed
340 341 342 343
{
    insideOut = io;
}

hrchilds's avatar
hrchilds committed
344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363
// ****************************************************************************
//  Method:  vtkVisItClipper::SetRemoveWholeCells
//
//  Purpose:
//    Tell the clipper if you want it to treat cells as atomic, and
//    simply remove any cell not entirely within the region.
//
//  Arguments:
//    lcw        the new setting
//
//  Programmer:  Jeremy Meredith
//  Creation:    August 29, 2006
//
// ****************************************************************************
void
vtkVisItClipper::SetRemoveWholeCells(bool rwc)
{
    removeWholeCells = rwc;
}

hrchilds's avatar
hrchilds committed
364 365 366 367 368
vtkUnstructuredGrid*
vtkVisItClipper::GetOtherOutput()
{
    return otherOutput;
}
hrchilds's avatar
hrchilds committed
369

hrchilds's avatar
hrchilds committed
370
// ****************************************************************************
hrchilds's avatar
hrchilds committed
371
//  Method:  vtkVisItClipper::Execute
hrchilds's avatar
hrchilds committed
372 373
//
//  Purpose:
374
//    Main execution method.  
hrchilds's avatar
hrchilds committed
375 376 377 378 379
//
//  Arguments:
//    none
//
//  Programmer:  Jeremy Meredith
380
//  Creation:    February 24, 2010
hrchilds's avatar
hrchilds committed
381 382
//
//  Modifications:
383 384 385
//    Jeremy Meredith, Wed Feb 24 10:18:33 EST 2010
//    Initial creation: unified the old rectilinear, structured, unstructured,
//    and polydata execution functions into this single function.
hrchilds's avatar
hrchilds committed
386 387
//
// ****************************************************************************
hrchilds's avatar
hrchilds committed
388
void
hrchilds's avatar
hrchilds committed
389
vtkVisItClipper::Execute()
hrchilds's avatar
hrchilds committed
390
{
391
    vtkDataSet *ds = GetInput();
hrchilds's avatar
hrchilds committed
392

393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426
    int do_type = ds->GetDataObjectType();
    vtkRectilinearGrid   *rg = NULL;
    vtkStructuredGrid    *sg = NULL;
    vtkUnstructuredGrid  *ug = NULL;
    vtkPolyData          *pg = NULL;

    // coordinate arrays for any mesh type
    float      *X       = NULL;
    float      *Y       = NULL;
    float      *Z       = NULL;
    float      *pts_ptr = NULL;

    // dimensions for structured grids
    int pt_dims[3];
    int cell_dims[3];
    int strideY;
    int strideZ;
    int ptstrideY;
    int ptstrideZ;

    // indices to convert structured grid cells to hexahedron/quadrilateral
    const int X_val[8] = { 0, 1, 1, 0, 0, 1, 1, 0 };
    const int Y_val[8] = { 0, 0, 1, 1, 0, 0, 1, 1 };
    const int Z_val[8] = { 0, 0, 0, 0, 1, 1, 1, 1 };

    // Set general input/output data
    int                  nCells = ds->GetNumberOfCells();
    vtkCellData         *inCD   = ds->GetCellData();
    vtkPointData        *inPD   = ds->GetPointData();
    vtkUnstructuredGrid *output = (vtkUnstructuredGrid*)GetOutput();
    vtkUnstructuredGrid *stuff_I_cant_clip = vtkUnstructuredGrid::New();

    bool twoD = false;
    if (do_type == VTK_RECTILINEAR_GRID || do_type == VTK_STRUCTURED_GRID)
hrchilds's avatar
hrchilds committed
427
    {
428 429 430 431 432 433 434 435 436 437 438 439 440 441 442
        if (do_type == VTK_RECTILINEAR_GRID)
        {
            rg = (vtkRectilinearGrid*)ds;
            rg->GetDimensions(pt_dims);
            X = (float* ) rg->GetXCoordinates()->GetVoidPointer(0);
            Y = (float* ) rg->GetYCoordinates()->GetVoidPointer(0);
            Z = (float* ) rg->GetZCoordinates()->GetVoidPointer(0);
        }
        else // do_type == VTK_STRUCTURED_GRID
        {
            sg = (vtkStructuredGrid*)ds;
            sg->GetDimensions(pt_dims);
            pts_ptr = (float*)sg->GetPoints()->GetVoidPointer(0);
        }

443
        twoD = (pt_dims[2] <= 1);
444 445 446 447 448 449 450
        cell_dims[0] = pt_dims[0]-1;
        cell_dims[1] = pt_dims[1]-1;
        cell_dims[2] = pt_dims[2]-1;
        strideY = cell_dims[0];
        strideZ = cell_dims[0]*cell_dims[1];
        ptstrideY = pt_dims[0];
        ptstrideZ = pt_dims[0]*pt_dims[1];
hrchilds's avatar
hrchilds committed
451 452 453
    }
    else if (do_type == VTK_UNSTRUCTURED_GRID)
    {
454 455 456 457 458
        ug = (vtkUnstructuredGrid*)ds;
        pts_ptr = (float*)ug->GetPoints()->GetVoidPointer(0);
        stuff_I_cant_clip->SetPoints(ug->GetPoints());
        stuff_I_cant_clip->GetPointData()->ShallowCopy(ug->GetPointData());
        stuff_I_cant_clip->Allocate(nCells);
hrchilds's avatar
hrchilds committed
459
    }
hrchilds's avatar
hrchilds committed
460 461
    else if (do_type == VTK_POLY_DATA)
    {
462 463 464 465 466 467
        pg = (vtkPolyData*)ds;
        pts_ptr = (float*)pg->GetPoints()->GetVoidPointer(0);
        stuff_I_cant_clip->SetPoints(pg->GetPoints());
        stuff_I_cant_clip->GetPointData()->ShallowCopy(pg->GetPointData());
        stuff_I_cant_clip->Allocate(nCells);

hrchilds's avatar
hrchilds committed
468
    }
hrchilds's avatar
hrchilds committed
469 470
    else
    {
471 472
        debug1 << "vtkVisItClipper: Can't operate on this dataset,\n";
        debug1 << "                 reverting to raw VTK code.\n";
hrchilds's avatar
hrchilds committed
473 474 475 476
        GeneralExecute();
    }

    int ptSizeGuess = (CellList == NULL
hrchilds's avatar
hrchilds committed
477 478
                         ? (int) pow(float(nCells), 0.6667f) * 5 + 100
                         : CellListSize*5 + 100);
hrchilds's avatar
hrchilds committed
479

480 481
    vtkVolumeFromVolume vfvIn(ds->GetNumberOfPoints(), ptSizeGuess);
    vtkVolumeFromVolume vfvOut(ds->GetNumberOfPoints(), ptSizeGuess);
482
    vtkVolumeFromVolume *useVFV;
hrchilds's avatar
hrchilds committed
483

484 485 486 487 488
    const int max_pts = 8;
    int cellType = twoD ? VTK_QUAD : VTK_HEXAHEDRON; // constant for struct grd
    int nCellPts = twoD ? 4 : 8;                     // constant for struct grd
    vtkIdType cellPtsStruct[8];
    vtkIdType *cellPts = cellPtsStruct; // for struct grd, we'll fill it
hrchilds's avatar
hrchilds committed
489 490

    int nToProcess = (CellList != NULL ? CellListSize : nCells);
491 492
    int numIcantClip = 0;
    for (int i = 0 ; i < nToProcess ; i++)
hrchilds's avatar
hrchilds committed
493
    {
494
        // Get the cell details
hrchilds's avatar
hrchilds committed
495
        int cellId = (CellList != NULL ? CellList[i] : i);
496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547
        int cellI = -1;
        int cellJ = -1;
        int cellK = -1;
        if (ug)
        {
            cellType = ug->GetCellType(cellId);
            ug->GetCellPoints(cellId, nCellPts, cellPts);
            // don't need cellI/J/K
        }
        else if (pg)
        {
            cellType = pg->GetCellType(cellId);
            pg->GetCellPoints(cellId, nCellPts, cellPts);
            // don't need cellI/J/K
        }
        else // structured grid
        {
            // cellType already set
            // nCellPts already set
            cellI = cellId % cell_dims[0];
            cellJ = (cellId/strideY) % cell_dims[1];
            cellK = (cellId/strideZ);
            for (int j = 0; j<nCellPts; j++)
            {
                cellPts[j] = (cellI + X_val[j]) +
                             (cellJ + Y_val[j])*ptstrideY +
                             (cellK + Z_val[j])*ptstrideZ;
            }
        }

        // If it's something we can't clip, save it for later
        bool canClip = false;
        switch (cellType)
        {
          case VTK_TETRA:
          case VTK_PYRAMID:
          case VTK_WEDGE:
          case VTK_HEXAHEDRON:
          case VTK_VOXEL:
          case VTK_TRIANGLE:
          case VTK_QUAD:
          case VTK_PIXEL:
          case VTK_LINE:
          case VTK_VERTEX:
            canClip = true;
            break;

          default:
            canClip = false;
            break;
        }
        if (!canClip)
hrchilds's avatar
hrchilds committed
548
        {
549 550 551 552 553 554 555 556 557 558
            if (numIcantClip == 0)
                stuff_I_cant_clip->GetCellData()->
                                       CopyAllocate(ds->GetCellData(), nCells);

            stuff_I_cant_clip->InsertNextCell(cellType, nCellPts, cellPts);
            stuff_I_cant_clip->GetCellData()->
                            CopyData(ds->GetCellData(), cellId, numIcantClip);
            numIcantClip++;
            continue;
        }
hrchilds's avatar
hrchilds committed
559

560 561 562 563 564
        // fill the dist functions and calculate lookup case
        int lookup_case = 0;
        float dist[max_pts];
        for (int j = nCellPts-1 ; j >= 0 ; j--)
        {
hrchilds's avatar
hrchilds committed
565 566
            if (clipFunction)
            {
567 568 569 570 571 572 573 574 575 576 577 578
                float ptRect[3];
                float *pt = ptRect;
                if (pts_ptr)
                {
                    pt = pts_ptr + 3*cellPts[j];
                }
                else
                {
                    pt[0] = X[cellI + X_val[j]];
                    pt[1] = Y[cellJ + Y_val[j]];
                    pt[2] = Z[cellK + Z_val[j]];
                }
hrchilds's avatar
hrchilds committed
579 580 581 582
                dist[j] = clipFunction->EvaluateFunction(pt[0],pt[1],pt[2]);
            }
            else // if (scalarArray)
            {
583
                float val = scalarArray[cellPts[j]];
hrchilds's avatar
hrchilds committed
584
                dist[j] = scalarCutoff - val;
hrchilds's avatar
hrchilds committed
585 586
            }

hrchilds's avatar
hrchilds committed
587 588 589 590 591
            if (dist[j] >= 0)
                lookup_case++;
            if (j > 0)
                lookup_case *= 2;
        }
hrchilds's avatar
hrchilds committed
592 593 594 595

        if (removeWholeCells && lookup_case != 0)
            lookup_case = ((1 << nCellPts) - 1);

596 597 598 599 600 601 602
        unsigned char  *splitCase = NULL;
        int             numOutput = 0;
        typedef int     edgeIndices[2];
        edgeIndices    *vertices_from_edges = NULL;

        int startIndex;
        switch (cellType)
hrchilds's avatar
hrchilds committed
603
        {
604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663
          case VTK_TETRA:
            startIndex = startClipShapesTet[lookup_case];
            splitCase  = &clipShapesTet[startIndex];
            numOutput  = numClipShapesTet[lookup_case];
            vertices_from_edges = tetVerticesFromEdges;
            break;
          case VTK_PYRAMID:
            startIndex = startClipShapesPyr[lookup_case];
            splitCase  = &clipShapesPyr[startIndex];
            numOutput  = numClipShapesPyr[lookup_case];
            vertices_from_edges = pyramidVerticesFromEdges;
            break;
          case VTK_WEDGE:
            startIndex = startClipShapesWdg[lookup_case];
            splitCase  = &clipShapesWdg[startIndex];
            numOutput  = numClipShapesWdg[lookup_case];
            vertices_from_edges = wedgeVerticesFromEdges;
            break;
          case VTK_HEXAHEDRON:
            startIndex = startClipShapesHex[lookup_case];
            splitCase  = &clipShapesHex[startIndex];
            numOutput  = numClipShapesHex[lookup_case];
            vertices_from_edges = hexVerticesFromEdges;
            break;
          case VTK_VOXEL:
            startIndex = startClipShapesVox[lookup_case];
            splitCase  = &clipShapesVox[startIndex];
            numOutput  = numClipShapesVox[lookup_case];
            vertices_from_edges = voxVerticesFromEdges;
            break;
          case VTK_TRIANGLE:
            startIndex = startClipShapesTri[lookup_case];
            splitCase  = &clipShapesTri[startIndex];
            numOutput  = numClipShapesTri[lookup_case];
            vertices_from_edges = triVerticesFromEdges;
            break;
          case VTK_QUAD:
            startIndex = startClipShapesQua[lookup_case];
            splitCase  = &clipShapesQua[startIndex];
            numOutput  = numClipShapesQua[lookup_case];
            vertices_from_edges = quadVerticesFromEdges;
            break;
          case VTK_PIXEL:
            startIndex = startClipShapesPix[lookup_case];
            splitCase  = &clipShapesPix[startIndex];
            numOutput  = numClipShapesPix[lookup_case];
            vertices_from_edges = pixelVerticesFromEdges;
            break;
          case VTK_LINE:
            startIndex = startClipShapesLin[lookup_case];
            splitCase  = &clipShapesLin[startIndex];
            numOutput  = numClipShapesLin[lookup_case];
            vertices_from_edges = lineVerticesFromEdges;
            break;
          case VTK_VERTEX:
            startIndex = startClipShapesVtx[lookup_case];
            splitCase  = &clipShapesVtx[startIndex];
            numOutput  = numClipShapesVtx[lookup_case];
            vertices_from_edges = NULL;
            break;
hrchilds's avatar
hrchilds committed
664 665
        }

666 667
        int            interpIDsIn[4];
        int            interpIDsOut[4];
668
        for (int j = 0 ; j < numOutput ; j++)
hrchilds's avatar
hrchilds committed
669 670 671 672
        {
            unsigned char shapeType = *splitCase++;
            {
                int npts;
673
                int interpID = -1;
hrchilds's avatar
hrchilds committed
674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692
                int color    = -1;
                switch (shapeType)
                {
                  case ST_HEX:
                    npts = 8;
                    color = *splitCase++;
                    break;
                  case ST_WDG:
                    npts = 6;
                    color = *splitCase++;
                    break;
                  case ST_PYR:
                    npts = 5;
                    color = *splitCase++;
                    break;
                  case ST_TET:
                    npts = 4;
                    color = *splitCase++;
                    break;
hrchilds's avatar
hrchilds committed
693 694 695 696 697 698 699 700
                  case ST_QUA:
                    npts = 4;
                    color = *splitCase++;
                    break;
                  case ST_TRI:
                    npts = 3;
                    color = *splitCase++;
                    break;
hrchilds's avatar
hrchilds committed
701 702 703 704 705 706 707 708
                  case ST_LIN:
                    npts = 2;
                    color = *splitCase++;
                    break;
                  case ST_VTX:
                    npts = 1;
                    color = *splitCase++;
                    break;
hrchilds's avatar
hrchilds committed
709
                  case ST_PNT:
710
                    interpID = *splitCase++;
hrchilds's avatar
hrchilds committed
711
                    color    = *splitCase++;
hrchilds's avatar
hrchilds committed
712 713 714 715 716
                    npts     = *splitCase++;
                    break;
                  default:
                    EXCEPTION1(ImproperUseException,
                               "An invalid output shape was found in "
hrchilds's avatar
hrchilds committed
717
                               "the ClipCases.");
hrchilds's avatar
hrchilds committed
718 719
                }

720
                useVFV = &vfvIn;
hrchilds's avatar
hrchilds committed
721 722 723
                if ((!insideOut && color == COLOR0) ||
                    ( insideOut && color == COLOR1))
                {
724 725 726 727 728 729 730 731 732 733
                    if (computeInsideAndOut)
                    {
                        useVFV = &vfvOut;
                    }
                    else
                    {
                        // We don't want this one; it's the wrong side.
                        splitCase += npts;
                        continue;
                    }
hrchilds's avatar
hrchilds committed
734 735 736 737 738 739
                }

                int shape[8];
                for (int p = 0 ; p < npts ; p++)
                {
                    unsigned char pt = *splitCase++;
hrchilds's avatar
hrchilds committed
740
                    if (pt <= P7)
hrchilds's avatar
hrchilds committed
741
                    {
hrchilds's avatar
hrchilds committed
742 743 744
                        // We know pt P0 must be >P0 since we already
                        // assume P0 == 0.  This is why we do not
                        // bother subtracting P0 from pt here.
745
                        shape[p] = cellPts[pt];
hrchilds's avatar
hrchilds committed
746 747 748
                    }
                    else if (pt >= EA && pt <= EL)
                    {
749 750
                        int pt1 = vertices_from_edges[pt-EA][0];
                        int pt2 = vertices_from_edges[pt-EA][1];
hrchilds's avatar
hrchilds committed
751 752 753 754 755 756 757 758 759
                        if (pt2 < pt1)
                        {
                            int tmp = pt2;
                            pt2 = pt1;
                            pt1 = tmp;
                        }
                        float dir = dist[pt2] - dist[pt1];
                        float amt = 0. - dist[pt1];
                        float percent = 1. - (amt / dir);
hrchilds's avatar
hrchilds committed
760 761 762 763

                        // We may have physically (though not logically)
                        // degenerate cells if percent==0 or percent==1.
                        // We could pretty easily and mostly safely clamp
764 765 766
                        // percent to the range [1e-4, 1. - 1e-4] here.
                        int ptId1 = cellPts[pt1];
                        int ptId2 = cellPts[pt2];
hrchilds's avatar
hrchilds committed
767 768 769

                        // deal with exact zero crossings if requested
                        if (clipFunction && useZeroCrossings)
770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793
                        {
                            if (pts_ptr)
                            {
                                AdjustPercentToZeroCrossing(pts_ptr,
                                                            ptId1, ptId2,
                                                            clipFunction,
                                                            &percent);
                            }
                            else
                            {
                                // fake a little points array for rgrids
                                float pt[6];
                                pt[0] = X[cellI + X_val[pt1]];
                                pt[1] = Y[cellJ + Y_val[pt1]];
                                pt[2] = Z[cellK + Z_val[pt1]];
                                pt[3] = X[cellI + X_val[pt2]];
                                pt[4] = Y[cellJ + Y_val[pt2]];
                                pt[5] = Z[cellK + Z_val[pt2]];
                                AdjustPercentToZeroCrossing(pt,
                                                            0, 1,
                                                            clipFunction,
                                                            &percent);
                            }
                        }
hrchilds's avatar
hrchilds committed
794
                                
795
                        shape[p] = useVFV->AddPoint(ptId1, ptId2, percent);
hrchilds's avatar
hrchilds committed
796 797 798
                    }
                    else if (pt >= N0 && pt <= N3)
                    {
799 800 801 802
                        if (useVFV == &vfvIn)
                            shape[p] = interpIDsIn[pt - N0];
                        else
                            shape[p] = interpIDsOut[pt - N0];
hrchilds's avatar
hrchilds committed
803 804 805 806 807
                    }
                    else
                    {
                        EXCEPTION1(ImproperUseException,
                                   "An invalid output point value "
hrchilds's avatar
hrchilds committed
808
                                   "was found in the ClipCases.");
hrchilds's avatar
hrchilds committed
809 810 811 812 813 814
                    }
                }

                switch (shapeType)
                {
                  case ST_HEX:
815
                    useVFV->AddHex(cellId,
816 817
                                   shape[0], shape[1], shape[2], shape[3],
                                   shape[4], shape[5], shape[6], shape[7]);
hrchilds's avatar
hrchilds committed
818 819
                    break;
                  case ST_WDG:
820
                    useVFV->AddWedge(cellId,
821 822
                                     shape[0], shape[1], shape[2],
                                     shape[3], shape[4], shape[5]);
hrchilds's avatar
hrchilds committed
823 824
                    break;
                  case ST_PYR:
825
                    useVFV->AddPyramid(cellId, shape[0], shape[1],
826
                                       shape[2], shape[3], shape[4]);
hrchilds's avatar
hrchilds committed
827 828
                    break;
                  case ST_TET:
829
                    useVFV->AddTet(cellId, shape[0], shape[1], shape[2], shape[3]);
hrchilds's avatar
hrchilds committed
830
                    break;
hrchilds's avatar
hrchilds committed
831
                  case ST_QUA:
832
                    useVFV->AddQuad(cellId, shape[0], shape[1], shape[2], shape[3]);
hrchilds's avatar
hrchilds committed
833 834
                    break;
                  case ST_TRI:
835
                    useVFV->AddTri(cellId, shape[0], shape[1], shape[2]);
hrchilds's avatar
hrchilds committed
836
                    break;
hrchilds's avatar
hrchilds committed
837
                  case ST_LIN:
838
                    useVFV->AddLine(cellId, shape[0], shape[1]);
hrchilds's avatar
hrchilds committed
839 840
                    break;
                  case ST_VTX:
841
                    useVFV->AddVertex(cellId, shape[0]);
hrchilds's avatar
hrchilds committed
842
                    break;
hrchilds's avatar
hrchilds committed
843
                  case ST_PNT:
844 845 846
                    interpIDsIn[interpID] = vfvIn.AddCentroidPoint(npts, shape);
                    if (computeInsideAndOut)
                        interpIDsOut[interpID] = vfvOut.AddCentroidPoint(npts, shape);
hrchilds's avatar
hrchilds committed
847 848 849 850 851 852
                    break;
                }
            }
        }
    }

853
    if (numIcantClip > 0)
854
    {
855 856 857 858 859 860 861 862
        vtkUnstructuredGrid *not_from_zoo  = vtkUnstructuredGrid::New();
        ClipDataset(stuff_I_cant_clip, not_from_zoo);
        
        vtkUnstructuredGrid *just_from_zoo = vtkUnstructuredGrid::New();
        if (pts_ptr)
            vfvIn.ConstructDataSet(inPD, inCD, just_from_zoo, pts_ptr);
        else
            vfvIn.ConstructDataSet(inPD, inCD, just_from_zoo, pt_dims,X,Y,Z);
hrchilds's avatar
hrchilds committed
863

864 865 866 867
        vtkAppendFilter *appender = vtkAppendFilter::New();
        appender->AddInput(not_from_zoo);
        appender->AddInput(just_from_zoo);
        appender->GetOutput()->Update();
hrchilds's avatar
hrchilds committed
868

869
        output->ShallowCopy(appender->GetOutput());
hrchilds's avatar
hrchilds committed
870

871
        if (computeInsideAndOut)
hrchilds's avatar
hrchilds committed
872
        {
873 874
            appender->RemoveInput(just_from_zoo);
            just_from_zoo->Delete();
hrchilds's avatar
hrchilds committed
875

876 877 878 879 880
            just_from_zoo = vtkUnstructuredGrid::New();
            if (pts_ptr)
                vfvOut.ConstructDataSet(inPD, inCD, just_from_zoo, pts_ptr);
            else
                vfvOut.ConstructDataSet(inPD, inCD, just_from_zoo, pt_dims,X,Y,Z);
hrchilds's avatar
hrchilds committed
881

882 883
            appender->AddInput(just_from_zoo);
            appender->GetOutput()->Update();
hrchilds's avatar
hrchilds committed
884

885 886 887
            if (otherOutput) otherOutput->Delete();
            otherOutput = vtkUnstructuredGrid::New();
            otherOutput->ShallowCopy(appender->GetOutput());
hrchilds's avatar
hrchilds committed
888
        }
889 890 891 892 893 894 895 896 897

        appender->Delete();
        just_from_zoo->Delete();
        not_from_zoo->Delete();
    }
    else
    {
        if (pts_ptr)
            vfvIn.ConstructDataSet(inPD, inCD, output, pts_ptr);
hrchilds's avatar
hrchilds committed
898
        else
899 900 901
            vfvIn.ConstructDataSet(inPD, inCD, output, pt_dims,X,Y,Z);

        if (computeInsideAndOut)
hrchilds's avatar
hrchilds committed
902
        {
903 904 905 906 907 908
            if (otherOutput) otherOutput->Delete();
            otherOutput = vtkUnstructuredGrid::New();
            if (pts_ptr)
                vfvOut.ConstructDataSet(inPD, inCD, otherOutput, pts_ptr);
            else
                vfvOut.ConstructDataSet(inPD, inCD, otherOutput, pt_dims,X,Y,Z);
hrchilds's avatar
hrchilds committed
909
        }
910
    }
hrchilds's avatar
hrchilds committed
911 912 913 914

    stuff_I_cant_clip->Delete();
}

hrchilds's avatar
hrchilds committed
915

hrchilds's avatar
hrchilds committed
916
void vtkVisItClipper::PrintSelf(ostream& os, vtkIndent indent)
hrchilds's avatar
hrchilds committed
917 918 919 920
{
    Superclass::PrintSelf(os,indent);
}

hrchilds's avatar
hrchilds committed
921
void vtkVisItClipper::GeneralExecute(void)
hrchilds's avatar
hrchilds committed
922 923 924 925
{
    ClipDataset(GetInput(), (vtkUnstructuredGrid*)GetOutput());
}

hrchilds's avatar
hrchilds committed
926 927 928 929 930 931 932
// ****************************************************************************
//  Modifications:
//
//    Hank Childs, Sat Mar 27 10:56:08 PST 2004
//    Work-around some funniness with VTK memory management.  (the funniness
//    is a bug with the vtkClipDataSet filter.)
//
933 934 935
//    Hank Childs, Sat Oct  6 15:37:11 PDT 2007
//    Fix bug with setting "inverse" for isovoluming.
//
hrchilds's avatar
hrchilds committed
936 937
// ****************************************************************************

hrchilds's avatar
hrchilds committed
938
void vtkVisItClipper::ClipDataset(vtkDataSet *in_ds,
939
                                  vtkUnstructuredGrid *out_ds)
hrchilds's avatar
hrchilds committed
940 941 942
{
    vtkClipDataSet *clipData = vtkClipDataSet::New();
    clipData->SetInput(in_ds);
943 944 945 946
    if (clipFunction)
    {
        clipData->SetClipFunction(clipFunction);
        clipData->GenerateClipScalarsOff();
947
        clipData->SetInsideOut(insideOut);
948 949 950 951 952 953 954
    }
    else
    {
        clipData->SetClipFunction(NULL);
        in_ds->GetPointData()->SetScalars(scalarArrayAsVTK);
        clipData->GenerateClipScalarsOff();
        clipData->SetValue(scalarCutoff);
955
        clipData->SetInsideOut(!insideOut);
956
    }
hrchilds's avatar
hrchilds committed
957
    clipData->Update();
hrchilds's avatar
hrchilds committed
958
    out_ds->ShallowCopy(clipData->GetOutput());
hrchilds's avatar
hrchilds committed
959 960 961
    clipData->Delete();
}