vtkVisItClipper.C 33.2 KB
Newer Older
hrchilds's avatar
hrchilds committed
1 2
/*****************************************************************************
*
3
* Copyright (c) 2000 - 2010, 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 388 389
//    Jeremy Meredith, Thu Feb 25 11:08:03 EST 2010
//    Don't forget to exit early if we have a dataset we can't understand.
//
390 391 392 393 394
//    Jeremy Meredith, Thu Feb 25 15:14:28 EST 2010
//    Allowing clipFunction usage to precalculate (most) values.  This
//    saves a good chunk of time in this mode since we were re-calculating
//    these values a number of times.
//
hrchilds's avatar
hrchilds committed
395
// ****************************************************************************
hrchilds's avatar
hrchilds committed
396
void
hrchilds's avatar
hrchilds committed
397
vtkVisItClipper::Execute()
hrchilds's avatar
hrchilds committed
398
{
399
    vtkDataSet *ds = GetInput();
hrchilds's avatar
hrchilds committed
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 427 428 429 430 431 432 433 434
    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
435
    {
436 437 438 439 440 441 442 443 444 445 446 447 448 449 450
        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);
        }

451
        twoD = (pt_dims[2] <= 1);
452 453 454 455 456 457 458
        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];
459 460
        // we can clip all of structured grids; don't set up
        // the "stuff_I_cant_clip" mesh
hrchilds's avatar
hrchilds committed
461 462 463
    }
    else if (do_type == VTK_UNSTRUCTURED_GRID)
    {
464 465 466 467 468
        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
469
    }
hrchilds's avatar
hrchilds committed
470 471
    else if (do_type == VTK_POLY_DATA)
    {
472 473 474 475 476 477
        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
478
    }
hrchilds's avatar
hrchilds committed
479 480
    else
    {
481 482
        debug1 << "vtkVisItClipper: Can't operate on this dataset,\n";
        debug1 << "                 reverting to raw VTK code.\n";
hrchilds's avatar
hrchilds committed
483
        GeneralExecute();
484
        return;
hrchilds's avatar
hrchilds committed
485 486
    }

487 488 489 490 491 492 493 494 495 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
    //
    // If we have a clip function, evaluate it on the points here
    // (better than doing it multiple times when evaluating each
    // node from each cell).
    //
    if (clipFunction)
    {
        int npts = ds->GetNumberOfPoints();
        scalarArray = new float[npts];
        scalarCutoff = 0;
        if (pts_ptr)
        {
            for (int i=0; i<npts; i++)
            {
                float *pt =  pts_ptr + 3*i;
                scalarArray[i] = 
                   -clipFunction->EvaluateFunction(pt[0],pt[1],pt[2]);
            }
        }
        else
        {
            int ctr = 0;
            for (int k=0; k<pt_dims[2]; k++)
            {
                for (int j=0; j<pt_dims[1]; j++)
                {
                    for (int i=0; i<pt_dims[0]; i++)
                    {
                        scalarArray[ctr++] = 
                           -clipFunction->EvaluateFunction(X[i],Y[j],Z[k]);
                    }
                }
            }
        }
    }


    //
    // Do the actual clipping here
    //
hrchilds's avatar
hrchilds committed
527
    int ptSizeGuess = (CellList == NULL
hrchilds's avatar
hrchilds committed
528 529
                         ? (int) pow(float(nCells), 0.6667f) * 5 + 100
                         : CellListSize*5 + 100);
hrchilds's avatar
hrchilds committed
530

531 532
    vtkVolumeFromVolume vfvIn(ds->GetNumberOfPoints(), ptSizeGuess);
    vtkVolumeFromVolume vfvOut(ds->GetNumberOfPoints(), ptSizeGuess);
533
    vtkVolumeFromVolume *useVFV;
hrchilds's avatar
hrchilds committed
534

535 536 537 538 539
    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
540 541

    int nToProcess = (CellList != NULL ? CellListSize : nCells);
542 543
    int numIcantClip = 0;
    for (int i = 0 ; i < nToProcess ; i++)
hrchilds's avatar
hrchilds committed
544
    {
545
        // Get the cell details
hrchilds's avatar
hrchilds committed
546
        int cellId = (CellList != NULL ? CellList[i] : i);
547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598
        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
599
        {
600
            if (numIcantClip == NULL)
601 602 603 604 605 606 607 608 609
                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
610

611 612 613 614 615
        // fill the dist functions and calculate lookup case
        int lookup_case = 0;
        float dist[max_pts];
        for (int j = nCellPts-1 ; j >= 0 ; j--)
        {
616 617
            float val = scalarArray[cellPts[j]];
            dist[j] = scalarCutoff - val;
hrchilds's avatar
hrchilds committed
618

hrchilds's avatar
hrchilds committed
619 620 621 622 623
            if (dist[j] >= 0)
                lookup_case++;
            if (j > 0)
                lookup_case *= 2;
        }
hrchilds's avatar
hrchilds committed
624 625 626 627

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

628 629 630 631 632 633 634
        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
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 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695
          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
696 697
        }

698 699
        int            interpIDsIn[4];
        int            interpIDsOut[4];
700
        for (int j = 0 ; j < numOutput ; j++)
hrchilds's avatar
hrchilds committed
701 702 703 704
        {
            unsigned char shapeType = *splitCase++;
            {
                int npts;
705
                int interpID = -1;
hrchilds's avatar
hrchilds committed
706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724
                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
725 726 727 728 729 730 731 732
                  case ST_QUA:
                    npts = 4;
                    color = *splitCase++;
                    break;
                  case ST_TRI:
                    npts = 3;
                    color = *splitCase++;
                    break;
hrchilds's avatar
hrchilds committed
733 734 735 736 737 738 739 740
                  case ST_LIN:
                    npts = 2;
                    color = *splitCase++;
                    break;
                  case ST_VTX:
                    npts = 1;
                    color = *splitCase++;
                    break;
hrchilds's avatar
hrchilds committed
741
                  case ST_PNT:
742
                    interpID = *splitCase++;
hrchilds's avatar
hrchilds committed
743
                    color    = *splitCase++;
hrchilds's avatar
hrchilds committed
744 745 746 747 748
                    npts     = *splitCase++;
                    break;
                  default:
                    EXCEPTION1(ImproperUseException,
                               "An invalid output shape was found in "
hrchilds's avatar
hrchilds committed
749
                               "the ClipCases.");
hrchilds's avatar
hrchilds committed
750 751
                }

752 753
                bool out = ((!insideOut && color == COLOR0) ||
                            ( insideOut && color == COLOR1));
754
                useVFV = &vfvIn;
755
                if (out)
hrchilds's avatar
hrchilds committed
756
                {
757 758 759 760 761 762 763 764 765 766
                    if (computeInsideAndOut)
                    {
                        useVFV = &vfvOut;
                    }
                    else
                    {
                        // We don't want this one; it's the wrong side.
                        splitCase += npts;
                        continue;
                    }
hrchilds's avatar
hrchilds committed
767 768 769 770 771 772
                }

                int shape[8];
                for (int p = 0 ; p < npts ; p++)
                {
                    unsigned char pt = *splitCase++;
hrchilds's avatar
hrchilds committed
773
                    if (pt <= P7)
hrchilds's avatar
hrchilds committed
774
                    {
hrchilds's avatar
hrchilds committed
775 776 777
                        // 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.
778
                        shape[p] = cellPts[pt];
hrchilds's avatar
hrchilds committed
779 780 781
                    }
                    else if (pt >= EA && pt <= EL)
                    {
782 783
                        int pt1 = vertices_from_edges[pt-EA][0];
                        int pt2 = vertices_from_edges[pt-EA][1];
hrchilds's avatar
hrchilds committed
784 785 786 787 788 789 790 791 792
                        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
793 794 795 796

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

                        // deal with exact zero crossings if requested
                        if (clipFunction && useZeroCrossings)
803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826
                        {
                            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
827
                                
828
                        shape[p] = useVFV->AddPoint(ptId1, ptId2, percent);
hrchilds's avatar
hrchilds committed
829 830 831
                    }
                    else if (pt >= N0 && pt <= N3)
                    {
832 833 834 835
                        if (useVFV == &vfvIn)
                            shape[p] = interpIDsIn[pt - N0];
                        else
                            shape[p] = interpIDsOut[pt - N0];
hrchilds's avatar
hrchilds committed
836 837 838 839 840
                    }
                    else
                    {
                        EXCEPTION1(ImproperUseException,
                                   "An invalid output point value "
hrchilds's avatar
hrchilds committed
841
                                   "was found in the ClipCases.");
hrchilds's avatar
hrchilds committed
842 843 844 845 846 847
                    }
                }

                switch (shapeType)
                {
                  case ST_HEX:
848
                    useVFV->AddHex(cellId,
849 850
                                   shape[0], shape[1], shape[2], shape[3],
                                   shape[4], shape[5], shape[6], shape[7]);
hrchilds's avatar
hrchilds committed
851 852
                    break;
                  case ST_WDG:
853
                    useVFV->AddWedge(cellId,
854 855
                                     shape[0], shape[1], shape[2],
                                     shape[3], shape[4], shape[5]);
hrchilds's avatar
hrchilds committed
856 857
                    break;
                  case ST_PYR:
858
                    useVFV->AddPyramid(cellId, shape[0], shape[1],
859
                                       shape[2], shape[3], shape[4]);
hrchilds's avatar
hrchilds committed
860 861
                    break;
                  case ST_TET:
862
                    useVFV->AddTet(cellId, shape[0], shape[1], shape[2], shape[3]);
hrchilds's avatar
hrchilds committed
863
                    break;
hrchilds's avatar
hrchilds committed
864
                  case ST_QUA:
865
                    useVFV->AddQuad(cellId, shape[0], shape[1], shape[2], shape[3]);
hrchilds's avatar
hrchilds committed
866 867
                    break;
                  case ST_TRI:
868
                    useVFV->AddTri(cellId, shape[0], shape[1], shape[2]);
hrchilds's avatar
hrchilds committed
869
                    break;
hrchilds's avatar
hrchilds committed
870
                  case ST_LIN:
871
                    useVFV->AddLine(cellId, shape[0], shape[1]);
hrchilds's avatar
hrchilds committed
872 873
                    break;
                  case ST_VTX:
874
                    useVFV->AddVertex(cellId, shape[0]);
hrchilds's avatar
hrchilds committed
875
                    break;
hrchilds's avatar
hrchilds committed
876
                  case ST_PNT:
877 878 879
                    interpIDsIn[interpID] = vfvIn.AddCentroidPoint(npts, shape);
                    if (computeInsideAndOut)
                        interpIDsOut[interpID] = vfvOut.AddCentroidPoint(npts, shape);
hrchilds's avatar
hrchilds committed
880 881 882 883 884 885
                    break;
                }
            }
        }
    }

886 887 888 889 890 891 892 893 894 895 896 897
    //
    // If we had a clip function, we created the scalarArray in this function
    //
    if (clipFunction)
    {
        delete[] scalarArray;
        scalarArray = NULL;
    }

    //
    // Construct the output data set.
    //
898
    if (numIcantClip > 0)
899
    {
900 901 902 903 904 905 906 907
        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
908

909 910 911 912
        vtkAppendFilter *appender = vtkAppendFilter::New();
        appender->AddInput(not_from_zoo);
        appender->AddInput(just_from_zoo);
        appender->GetOutput()->Update();
hrchilds's avatar
hrchilds committed
913

914
        output->ShallowCopy(appender->GetOutput());
hrchilds's avatar
hrchilds committed
915

916
        if (computeInsideAndOut)
hrchilds's avatar
hrchilds committed
917
        {
918 919
            appender->RemoveInput(just_from_zoo);
            just_from_zoo->Delete();
hrchilds's avatar
hrchilds committed
920

921 922 923 924 925
            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
926

927 928
            appender->AddInput(just_from_zoo);
            appender->GetOutput()->Update();
hrchilds's avatar
hrchilds committed
929

930 931 932
            if (otherOutput) otherOutput->Delete();
            otherOutput = vtkUnstructuredGrid::New();
            otherOutput->ShallowCopy(appender->GetOutput());
hrchilds's avatar
hrchilds committed
933
        }
934 935 936 937 938 939 940 941 942

        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
943
        else
944 945 946
            vfvIn.ConstructDataSet(inPD, inCD, output, pt_dims,X,Y,Z);

        if (computeInsideAndOut)
hrchilds's avatar
hrchilds committed
947
        {
948 949 950 951 952 953
            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
954
        }
955
    }
hrchilds's avatar
hrchilds committed
956 957 958 959

    stuff_I_cant_clip->Delete();
}

hrchilds's avatar
hrchilds committed
960

hrchilds's avatar
hrchilds committed
961
void vtkVisItClipper::PrintSelf(ostream& os, vtkIndent indent)
hrchilds's avatar
hrchilds committed
962 963 964 965
{
    Superclass::PrintSelf(os,indent);
}

hrchilds's avatar
hrchilds committed
966
void vtkVisItClipper::GeneralExecute(void)
hrchilds's avatar
hrchilds committed
967 968 969 970
{
    ClipDataset(GetInput(), (vtkUnstructuredGrid*)GetOutput());
}

hrchilds's avatar
hrchilds committed
971 972 973 974 975 976 977
// ****************************************************************************
//  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.)
//
978 979 980
//    Hank Childs, Sat Oct  6 15:37:11 PDT 2007
//    Fix bug with setting "inverse" for isovoluming.
//
hrchilds's avatar
hrchilds committed
981 982
// ****************************************************************************

hrchilds's avatar
hrchilds committed
983
void vtkVisItClipper::ClipDataset(vtkDataSet *in_ds,
984
                                  vtkUnstructuredGrid *out_ds)
hrchilds's avatar
hrchilds committed
985 986 987
{
    vtkClipDataSet *clipData = vtkClipDataSet::New();
    clipData->SetInput(in_ds);
988 989 990 991
    if (clipFunction)
    {
        clipData->SetClipFunction(clipFunction);
        clipData->GenerateClipScalarsOff();
992
        clipData->SetInsideOut(insideOut);
993 994 995 996 997 998 999
    }
    else
    {
        clipData->SetClipFunction(NULL);
        in_ds->GetPointData()->SetScalars(scalarArrayAsVTK);
        clipData->GenerateClipScalarsOff();
        clipData->SetValue(scalarCutoff);
1000
        clipData->SetInsideOut(!insideOut);
1001
    }
hrchilds's avatar
hrchilds committed
1002
    clipData->Update();
hrchilds's avatar
hrchilds committed
1003
    out_ds->ShallowCopy(clipData->GetOutput());
hrchilds's avatar
hrchilds committed
1004 1005 1006
    clipData->Delete();
}