[vtk-developers] vtkCutter modification to use vtkCellLocator

Rostislav Khlebnikov rostislav.khlebnikov at kcl.ac.uk
Thu Sep 18 15:14:29 EDT 2014


Hi,

the code is attached. As I said it's quite messy - just ask me if you 
can't figure out something. The functions I added start with small 
letters as opposed to the copy-pasted and modified from VTK which start 
with a capital letter.

Rostislav.


On 15/09/2014 19:28, Gerrick Bivins wrote:
>
> Yes. I would definitely be interested in seeing the code.
>
> Gerrick
>
> *From:*Rostislav Khlebnikov [mailto:rostislav.khlebnikov at kcl.ac.uk]
> *Sent:* Monday, September 15, 2014 12:07 PM
> *To:* Gerrick Bivins; Berk Geveci
> *Cc:* VTK Developers
> *Subject:* Re: [vtk-developers] vtkCutter modification to use 
> vtkCellLocator
>
> Hi,
>
> well, basically, I feel that it is necessary to create a proper change 
> to use Gerrit workflow successfully. Unfortunately, I don't quite have 
> the time to do this properly for several reasons. One of them is that 
> I actually don't need this change - for my project I integrated the 
> CGAL-based cutting approach with AABB_tree which outperforms my VTK 
> test - for example for the same test with 500k triangles CGAL-based 
> cutting takes only 2ms. This will soon be integrated to MITK (in case 
> anyone is interested, the pull request can be found here: 
> https://github.com/MITK/MITK/pull/73). Making a proper VTK change 
> would involve quite a lot of effort. Given that vtkCutter is very 
> general and rejecting the octree nodes might be complex for arbitrary 
> implicit functions, it seems to require quite a lot of work to avoid 
> incorrect usage.
>
> I can clean up the code a bit and post it to the mailing list if you 
> are interested.
>
> All best,
>    Rostislav.
>
> On 15/09/2014 15:02, Gerrick Bivins wrote:
>
>     Hi Rostislav,
>
>     Did anything happen with this?
>
>     I would definitely be interested in this
>     capability/feature/improvement.
>
>     Gerrick
>
>     *From:*vtk-developers [mailto:vtk-developers-bounces at vtk.org] *On
>     Behalf Of *Berk Geveci
>     *Sent:* Tuesday, August 12, 2014 11:44 AM
>     *To:* Rostislav Khlebnikov
>     *Cc:* VTK Developers
>     *Subject:* Re: [vtk-developers] vtkCutter modification to use
>     vtkCellLocator
>
>     Hi Rostislav,
>
>     This is great. Can you push the code to Gerrit
>     (http://review.source.kitware.com/)?
>
>     Best,
>
>     -berk
>
>     On Sat, Aug 9, 2014 at 2:52 PM, Rostislav Khlebnikov
>     <rostislav.khlebnikov at kcl.ac.uk
>     <mailto:rostislav.khlebnikov at kcl.ac.uk>> wrote:
>
>     Hi guys,
>
>     just thought that the VTK developers might be interested in the
>     test I made recently.
>
>     Bascially what I did is subclass the vtkCellLocator to implement
>     the visitor pattern. It allows to filter the octree nodes (e.g.
>     reject the octree nodes that are of no interest early) and for
>     leafs, to visit the cells of the polydata.
>     Then, I subclassed the vtkCutter to use this cell locator instead
>     of straight iteration over cells. Likely for very complex implicit
>     cutting functions this wouldnt be of any benefit. However, for
>     those which can test bboxes as being of interest or not, rejecting
>     whole subtrees, such as vtkPlane, it is very useful. I tested it
>     on a surface with 500k triangles and the cutter performance went
>     from 88ms per cut to 14ms.
>
>     The code I made is quite ugly due to limited extensibility of
>     vtkCutter, so I had to copy-paste the RequestData method and make
>     my changes instead of, say, overriding IterateCells method which
>     would just provide the cells to cut through.
>
>     But anyway, I was wondering if you guys are intersted in this kind
>     of functionality? Would you like me to send you the code?
>
>     All best,
>        Rostislav.
>
>     _______________________________________________
>     Powered by www.kitware.com <http://www.kitware.com>
>
>     Visit other Kitware open-source projects at
>     http://www.kitware.com/opensource/opensource.html
>     <http://www.kitware.com/opensource/opensource.html>
>
>     Follow this link to subscribe/unsubscribe:
>     http://public.kitware.com/mailman/listinfo/vtk-developers
>
>     ------------------------------------------------------------------------
>
>     This e-mail, including any attached files, may contain
>     confidential and privileged information for the sole use of the
>     intended recipient. Any review, use, distribution, or disclosure
>     by others is strictly prohibited. If you are not the intended
>     recipient (or authorized to receive information for the intended
>     recipient), please contact the sender by reply e-mail and delete
>     all copies of this message.
>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/vtk-developers/attachments/20140918/8075622e/attachment-0002.html>
-------------- next part --------------
#include <functional>
#include <queue>
#include <chrono>

#include <vtkSphereSource.h>

#include <vtkPolyData.h>
#include <vtkPolyDataReader.h>
#include <vtkCutter.h>
#include <vtkPlane.h>
#include <vtkIdList.h>
#include <vtkCell.h>

#include <vtkCellLocator.h>
#include <vtkCutter.h>
#include <vtkBox.h>
#include <vtkInformation.h>
#include <vtkInformationVector.h>
#include <vtkRectilinearGrid.h>
#include <vtkDoubleArray.h>
#include <vtkCellArray.h>
#include <vtkPointData.h>
#include <vtkCellData.h>
#include <vtkIncrementalPointLocator.h>
#include <vtkSmartPointer.h>
#include <vtkPolygonBuilder.h>
#include <vtkXMLPolyDataWriter.h>
#include <vtkCleanPolyData.h>
#include <vtkQuad.h>

class vtkContourHelper
{
public:
    vtkContourHelper(vtkIncrementalPointLocator *locator,
        vtkCellArray *verts,
        vtkCellArray *lines,
        vtkCellArray* polys,
        vtkPointData *inPd,
        vtkCellData *inCd,
        vtkPointData* outPd,
        vtkCellData *outCd,
        int estimatedSize,
        bool outputTriangles);
    ~vtkContourHelper();
    void Contour(vtkCell* cell, double value, vtkDataArray *cellScalars, vtkIdType cellId);

private:
    vtkIncrementalPointLocator* Locator;
    vtkCellArray* Verts;
    vtkCellArray* Lines;
    vtkCellArray* Polys;
    vtkPointData* InPd;
    vtkCellData* InCd;
    vtkPointData* OutPd;
    vtkCellData* OutCd;
    vtkSmartPointer<vtkCellData> TriOutCd;

    vtkCellArray* Tris;
    vtkPolygonBuilder PolyBuilder;
    vtkIdList* Poly;
    bool GenerateTriangles;
};

#include "vtkIncrementalPointLocator.h"
#include "vtkCellArray.h"
#include "vtkPointData.h"
#include "vtkCellData.h"
#include "vtkPolygonBuilder.h"
#include "vtkIdList.h"
#include "vtkCell.h"
#include "vtkDataArray.h"

vtkContourHelper::vtkContourHelper(vtkIncrementalPointLocator *locator,
    vtkCellArray *verts,
    vtkCellArray *lines,
    vtkCellArray* polys,
    vtkPointData *inPd,
    vtkCellData *inCd,
    vtkPointData* outPd,
    vtkCellData *outCd,
    int estimatedSize,
    bool outputTriangles) :
    Locator(locator),
    Verts(verts),
    Lines(lines),
    Polys(polys),
    InPd(inPd),
    InCd(inCd),
    OutPd(outPd),
    OutCd(outCd),
    GenerateTriangles(outputTriangles)
{
    this->Tris = vtkCellArray::New();
    this->TriOutCd = vtkCellData::New();
    if (this->GenerateTriangles)
    {
        this->Tris->Allocate(estimatedSize, estimatedSize / 2);
        this->TriOutCd->Initialize();
    }
    this->Poly = vtkIdList::New();
}

vtkContourHelper::~vtkContourHelper()
{
    this->Tris->Delete();
    this->TriOutCd->Delete();
    this->Poly->FastDelete();
}

void vtkContourHelper::Contour(vtkCell* cell, double value, vtkDataArray *cellScalars, vtkIdType cellId)
{
    bool mergeTriangles = (!this->GenerateTriangles) && cell->GetCellDimension() == 3;
    vtkCellData* outCD;
    vtkCellArray* outPoly;
    if (mergeTriangles)
    {
        outPoly = this->Tris;
        outCD = this->TriOutCd;
    }
    else
    {
        outPoly = this->Polys;
        outCD = this->OutCd;
    }
    cell->Contour(value, cellScalars, this->Locator, this->Verts, this->Lines,
        outPoly, this->InPd, this->OutPd, this->InCd, cellId, outCD);
    if (mergeTriangles)
    {
        this->PolyBuilder.Reset();

        vtkIdType cellSize;
        vtkIdType* cellVerts;
        while (this->Tris->GetNextCell(cellSize, cellVerts))
        {
            if (cellSize == 3)
            {
                this->PolyBuilder.InsertTriangle(cellVerts);
            }
            else //for whatever reason, the cell contouring is already outputing polys
            {
                vtkIdType outCellId = this->Polys->InsertNextCell(cellSize, cellVerts);
                this->OutCd->CopyData(this->InCd, cellId, outCellId);
            }
        }

        this->PolyBuilder.GetPolygon(this->Poly);
        if (this->Poly->GetNumberOfIds() != 0)
        {
            vtkIdType outCellId = this->Polys->InsertNextCell(this->Poly);
            this->OutCd->CopyData(this->InCd, cellId, outCellId);
        }
    }
}

class CellLocatorVisitor : public vtkCellLocator {
public:
    vtkTypeMacro(CellLocatorVisitor, vtkCellLocator);

    static CellLocatorVisitor *New();

    int getBoundsByIndex(int idx, double bounds[6])
    {
        int nSubdiv = 1;
        int level = 0;
        int nNodes = 1;

        while (idx >= nNodes) {
            nSubdiv *= 2;
            nNodes += nSubdiv * nSubdiv * nSubdiv;
            level++;
        }

        int firstInLevel = nNodes - nSubdiv*nSubdiv*nSubdiv;

        double sizeWithinLevel[3];
        for (int i = 0; i < 3; ++i) {
            sizeWithinLevel[i] = (this->Bounds[2 * i + 1] - this->Bounds[2 * i]) / nSubdiv;
        }

        idx -= firstInLevel;

        for (int i = 0; i < 3; ++i) {
            bounds[2 * i] = this->Bounds[2 * i] + sizeWithinLevel[i] * (idx % nSubdiv);
            bounds[2 * i + 1] = bounds[2 * i] + sizeWithinLevel[i];
            idx /= nSubdiv;
        }
        return level;
    }

    void getChildren(int idx, int children[8])
    {
        int nSubdiv = 1;
        int level = 0;
        int nNodes = 1;

        while (idx >= nNodes) {
            nSubdiv *= 2;
            nNodes += nSubdiv * nSubdiv * nSubdiv;
            level++;
        }

        int firstInParentLevel = nNodes - nSubdiv*nSubdiv*nSubdiv;

        idx -= firstInParentLevel;

        int firstChildIndex = nNodes;
        int nNodesInDirection = 2;
        for (int i = 0; i < 3; ++i) {
            firstChildIndex += (idx % (nSubdiv)) * nNodesInDirection;
            nNodesInDirection *= 2 * nSubdiv;
            idx /= (nSubdiv);
        }

        int outit = 0;
        for (int k = 0; k < 2; k++) {
            for (int j = 0; j < 2; j++) {
                for (int i = 0; i < 2; i++) {
                    children[outit++] = firstChildIndex + i + 2 * nSubdiv * j + 4 * nSubdiv * nSubdiv * k;
                }
            }
        }
    }
    void visitCells(std::function<bool(double[6])> cellChecker, std::function<void(vtkCell*, int)> visitor)
    {
        this->BuildLocatorIfNeeded();
        this->ClearCellHasBeenVisited();

        std::queue<int> octantsToVisit;

        octantsToVisit.push(0);

        while (!octantsToVisit.empty()) {
            int currentOctant = octantsToVisit.front();
            octantsToVisit.pop();

            // Get the octant bounds 
            double bounds[6];
            int level = getBoundsByIndex(currentOctant, bounds);

            if (cellChecker(bounds)) {
                vtkIdList* cellsIds = this->Tree[currentOctant];
                if (level == this->Level) {
                    // Leaf
                    if (!cellsIds) {
                        // Empty leaf
                        continue;
                    }

                    for (int i = 0; i < cellsIds->GetNumberOfIds(); ++i) {
                        int cellId = cellsIds->GetId(i);
                        if (!this->CellHasBeenVisited[cellId]) {
                            visitor(this->DataSet->GetCell(cellId), cellId);
                            this->CellHasBeenVisited[cellId] = 1;
                        }
                    }
                }
                else {
                    vtkIdList* cellsIds = this->Tree[currentOctant];
                    if (cellsIds != (void*)1) {
                        continue;
                    }

                    int children[8];

                    // Get children
                    getChildren(currentOctant, children);
                    for (int i = 0; i < 8; ++i) {
                        octantsToVisit.push(children[i]);
                    }
                }
            }
        }
    }
};

CellLocatorVisitor* CellLocatorVisitor::New()
{
    return new CellLocatorVisitor;
}

class CutterWithLocator : public vtkCutter {
public:
    vtkTypeMacro(CutterWithLocator, vtkCutter);

    static CutterWithLocator *New();

    void setLocator(CellLocatorVisitor* locator)
    {
        _locator = locator;
    }

    int RequestData(
        vtkInformation *request,
        vtkInformationVector **inputVector,
        vtkInformationVector *outputVector)
    {
        // get the info objects
        vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
        vtkInformation *outInfo = outputVector->GetInformationObject(0);

        // get the input and output
        vtkDataSet *input = vtkDataSet::SafeDownCast(
            inInfo->Get(vtkDataObject::DATA_OBJECT()));
        vtkPolyData *output = vtkPolyData::SafeDownCast(
            outInfo->Get(vtkDataObject::DATA_OBJECT()));

        vtkDebugMacro(<< "Executing cutter");
        if (!this->CutFunction)
        {
            vtkErrorMacro("No cut function specified");
            return 0;
        }

        if (!input)
        {
            // this could be a table in a multiblock structure, i.e. no cut!
            return 0;
        }

        if (input->GetNumberOfPoints() < 1 || this->GetNumberOfContours() < 1)
        {
            return 1;
        }

        if ((input->GetDataObjectType() == VTK_STRUCTURED_POINTS ||
            input->GetDataObjectType() == VTK_IMAGE_DATA) &&
            input->GetCell(0) && input->GetCell(0)->GetCellDimension() >= 3)
        {
            this->StructuredPointsCutter(input, output, request, inputVector, outputVector);
        }
        else if (input->GetDataObjectType() == VTK_STRUCTURED_GRID &&
            input->GetCell(0) &&
            input->GetCell(0)->GetCellDimension() >= 3)
        {
            this->StructuredGridCutter(input, output);
        }
        else if (input->GetDataObjectType() == VTK_RECTILINEAR_GRID &&
            static_cast<vtkRectilinearGrid *>(input)->GetDataDimension() == 3)
        {
            this->RectilinearGridCutter(input, output);
        }
        else if (input->GetDataObjectType() == VTK_UNSTRUCTURED_GRID)
        {
            vtkDebugMacro(<< "Executing Unstructured Grid Cutter");
            this->UnstructuredGridCutter(input, output);
        }
        else
        {
            vtkDebugMacro(<< "Executing DataSet Cutter");
            this->DataSetVisitorCutter(input, output);
        }

        return 1;
    }

    //----------------------------------------------------------------------------
    void DataSetVisitorCutter(vtkDataSet *input, vtkPolyData *output)
    {
        vtkIdType i;
        int iter;
        vtkPoints *cellPts;
        vtkDoubleArray *cellScalars;
        vtkCellArray *newVerts, *newLines, *newPolys;
        vtkPoints *newPoints;
        vtkDoubleArray *cutScalars;
        double value, s;
        vtkIdType estimatedSize, numCells = input->GetNumberOfCells();
        vtkIdType numPts = input->GetNumberOfPoints();
        int numCellPts;
        vtkPointData *inPD, *outPD;
        vtkCellData *inCD = input->GetCellData(), *outCD = output->GetCellData();
        vtkIdList *cellIds;
        int numContours = this->ContourValues->GetNumberOfContours();
        int abortExecute = 0;

        cellScalars = vtkDoubleArray::New();

        // Create objects to hold output of contour operation
        //
        estimatedSize = static_cast<vtkIdType>(
            pow(static_cast<double>(numCells), .75)) * numContours;
        estimatedSize = estimatedSize / 1024 * 1024; //multiple of 1024
        if (estimatedSize < 1024)
        {
            estimatedSize = 1024;
        }

        newPoints = vtkPoints::New();
        // set precision for the points in the output
        if (this->OutputPointsPrecision == vtkAlgorithm::DEFAULT_PRECISION)
        {
            vtkPointSet *inputPointSet = vtkPointSet::SafeDownCast(input);
            if (inputPointSet)
            {
                newPoints->SetDataType(inputPointSet->GetPoints()->GetDataType());
            }
            else
            {
                newPoints->SetDataType(VTK_FLOAT);
            }
        }
        else if (this->OutputPointsPrecision == vtkAlgorithm::SINGLE_PRECISION)
        {
            newPoints->SetDataType(VTK_FLOAT);
        }
        else if (this->OutputPointsPrecision == vtkAlgorithm::DOUBLE_PRECISION)
        {
            newPoints->SetDataType(VTK_DOUBLE);
        }
        newPoints->Allocate(estimatedSize, estimatedSize / 2);
        newVerts = vtkCellArray::New();
        newVerts->Allocate(estimatedSize, estimatedSize / 2);
        newLines = vtkCellArray::New();
        newLines->Allocate(estimatedSize, estimatedSize / 2);
        newPolys = vtkCellArray::New();
        newPolys->Allocate(estimatedSize, estimatedSize / 2);
        cutScalars = vtkDoubleArray::New();
        cutScalars->SetNumberOfTuples(numPts);

        // Interpolate data along edge. If generating cut scalars, do necessary setup
        if (this->GenerateCutScalars)
        {
            inPD = vtkPointData::New();
            inPD->ShallowCopy(input->GetPointData());//copies original attributes
            inPD->SetScalars(cutScalars);
        }
        else
        {
            inPD = input->GetPointData();
        }
        outPD = output->GetPointData();
        outPD->InterpolateAllocate(inPD, estimatedSize, estimatedSize / 2);
        outCD->CopyAllocate(inCD, estimatedSize, estimatedSize / 2);

        // locator used to merge potentially duplicate points
        if (this->Locator == NULL)
        {
            this->CreateDefaultLocator();
        }
        this->Locator->InitPointInsertion(newPoints, input->GetBounds());

        // Loop over all points evaluating scalar function at each point
        //
        for (i = 0; i < numPts; i++)
        {
            s = this->CutFunction->FunctionValue(input->GetPoint(i));
            cutScalars->SetComponent(i, 0, s);
        }

        // Compute some information for progress methods
        //
        vtkIdType numCuts = numContours*numCells;
        vtkIdType progressInterval = numCuts / 20 + 1;
        int cut = 0;

        vtkContourHelper helper(this->Locator, newVerts, newLines, newPolys, inPD, inCD, outPD, outCD, estimatedSize, this->GenerateTriangles != 0);
        if (true || this->SortBy == VTK_SORT_BY_CELL)
        {
            // Loop over all contour values.  Then for each contour value,
            // loop over all cells.
            //
            // This is going to have a problem if the input has 2D and 3D cells.
            // I am fixing a bug where cell data is scrambled becauses with
            // vtkPolyData output, verts and lines have lower cell ids than triangles.
            for (iter = 0; iter < numContours && !abortExecute; iter++)
            {
                _locator->visitCells(
                    [&](double bounds[6]) {
                    double sign = this->CutFunction->EvaluateFunction(bounds[0], bounds[2], bounds[4]);

                    bool oneSide = true;
                    for (int i = 0; i < 2 && oneSide; ++i) {
                        for (int j = 2; j < 4 && oneSide; ++j) {
                            for (int k = 4; k < 6; ++k) {
                                double sign2 = this->CutFunction->EvaluateFunction(bounds[i], bounds[j], bounds[k]);
                                if (sign * sign2 <= 0) {
                                    oneSide = false;
                                    break;
                                }
                            }
                        }
                    }

                    return !oneSide;
                },




                    [&](vtkCell* cell, int cellId) {
                    if (!(++cut % progressInterval))
                    {
                        vtkDebugMacro(<< "Cutting #" << cut);
                        this->UpdateProgress(static_cast<double>(cut) / numCuts);
                        abortExecute = this->GetAbortExecute();
                    }

                    cellPts = cell->GetPoints();
                    cellIds = cell->GetPointIds();

                    numCellPts = cellPts->GetNumberOfPoints();
                    cellScalars->SetNumberOfTuples(numCellPts);
                    for (i = 0; i < numCellPts; i++)
                    {
                        s = cutScalars->GetComponent(cellIds->GetId(i), 0);
                        cellScalars->SetTuple(i, &s);
                    }

                    value = this->ContourValues->GetValue(iter);

                    helper.Contour(cell, value, cellScalars, cellId);
                });
            } // for all contour values
        } // sort by cell

        //         else // VTK_SORT_BY_VALUE:
        //         {
        //             // Three passes over the cells to process lower dimensional cells first.
        //             // For poly data output cells need to be added in the order:
        //             // verts, lines and then polys, or cell data gets mixed up.
        //             // A better solution is to have an unstructured grid output.
        //             // I create a table that maps cell type to cell dimensionality,
        //             // because I need a fast way to get cell dimensionality.
        //             // This assumes GetCell is slow and GetCellType is fast.
        //             // I do not like hard coding a list of cell types here,
        //             // but I do not want to add GetCellDimension(vtkIdType cellId)
        //             // to the vtkDataSet API.  Since I anticipate that the output
        //             // will change to vtkUnstructuredGrid.  This temporary solution
        //             // is acceptable.
        //             //
        //             int cellType;
        //             unsigned char cellTypeDimensions[VTK_NUMBER_OF_CELL_TYPES];
        //             vtkCutter::GetCellTypeDimensions(cellTypeDimensions);
        //             int dimensionality;
        //             // We skip 0d cells (points), because they cannot be cut (generate no data).
        //             for (dimensionality = 1; dimensionality <= 3; ++dimensionality)
        //             {
        //                 // Loop over all cells; get scalar values for all cell points
        //                 // and process each cell.
        //                 //
        //                 for (cellId = 0; cellId < numCells && !abortExecute; cellId++)
        //                 {
        //                     // I assume that "GetCellType" is fast.
        //                     cellType = input->GetCellType(cellId);
        //                     if (cellType >= VTK_NUMBER_OF_CELL_TYPES)
        //                     { // Protect against new cell types added.
        //                         vtkErrorMacro("Unknown cell type " << cellType);
        //                         continue;
        //                     }
        //                     if (cellTypeDimensions[cellType] != dimensionality)
        //                     {
        //                         continue;
        //                     }
        //                     input->GetCell(cellId, cell);
        //                     cellPts = cell->GetPoints();
        //                     cellIds = cell->GetPointIds();
        // 
        //                     numCellPts = cellPts->GetNumberOfPoints();
        //                     cellScalars->SetNumberOfTuples(numCellPts);
        //                     for (i = 0; i < numCellPts; i++)
        //                     {
        //                         s = cutScalars->GetComponent(cellIds->GetId(i), 0);
        //                         cellScalars->SetTuple(i, &s);
        //                     }
        // 
        //                     // Loop over all contour values.
        //                     for (iter = 0; iter < numContours && !abortExecute; iter++)
        //                     {
        //                         if (dimensionality == 3 && !(++cut % progressInterval))
        //                         {
        //                             vtkDebugMacro(<< "Cutting #" << cut);
        //                             this->UpdateProgress(static_cast<double>(cut) / numCuts);
        //                             abortExecute = this->GetAbortExecute();
        //                         }
        //                         value = this->ContourValues->GetValue(iter);
        //                         helper.Contour(cell, value, cellScalars, cellId);
        //                     } // for all contour values
        //                 } // for all cells
        //             } // for all dimensions.
        //         } // sort by value

        // Update ourselves.  Because we don't know upfront how many verts, lines,
        // polys we've created, take care to reclaim memory.
        //
        cellScalars->Delete();
        cutScalars->Delete();

        if (this->GenerateCutScalars)
        {
            inPD->Delete();
        }

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

        if (newVerts->GetNumberOfCells())
        {
            output->SetVerts(newVerts);
        }
        newVerts->Delete();

        if (newLines->GetNumberOfCells())
        {
            output->SetLines(newLines);
        }
        newLines->Delete();

        if (newPolys->GetNumberOfCells())
        {
            output->SetPolys(newPolys);
        }
        newPolys->Delete();

        this->Locator->Initialize();//release any extra memory
        output->Squeeze();
    }

    CellLocatorVisitor* _locator;
};

CutterWithLocator* CutterWithLocator::New()
{
    return new CutterWithLocator;
}

int main()
{
    auto sphereSource = vtkSmartPointer<vtkSphereSource>::New();
    sphereSource->SetPhiResolution(300);
    sphereSource->SetThetaResolution(500);
    sphereSource->Update();

    auto plane = vtkSmartPointer<vtkPlane>::New();
    plane->SetOrigin(0.2, 0.1, -0.1);

    vtkCutter* cutter = vtkCutter::New();
    cutter->SetCutFunction(plane);
    cutter->SetInputData(sphereSource->GetOutput());
    cutter->GenerateCutScalarsOff();

    int nIter = 100;

    auto start = std::chrono::high_resolution_clock::now();

    for (int i = 0; i < nIter; ++i) {
        cutter->Modified();
        cutter->Update();
    }

    auto end = std::chrono::high_resolution_clock::now();
    auto dur = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);

    std::cout << "Normal vtkCutter " << dur.count() / (float)nIter << "ms, " << cutter->GetOutput()->GetNumberOfCells() << " cells\n";

    CellLocatorVisitor* locator = CellLocatorVisitor::New();
    locator->SetDataSet(sphereSource->GetOutput());
    locator->BuildLocator();

    CutterWithLocator* cutter2 = CutterWithLocator::New();
    cutter2->setLocator(locator);
    cutter2->SetCutFunction(plane);
    cutter2->SetInputData(sphereSource->GetOutput());
    cutter2->GenerateCutScalarsOff();

    start = std::chrono::high_resolution_clock::now();

    for (int i = 0; i < nIter; ++i) {
        cutter2->Modified();
        cutter2->Update();
    }

    end = std::chrono::high_resolution_clock::now();
    dur = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);

    std::cout << "'New' cutter " << dur.count() / (float)nIter << "ms, " << cutter2->GetOutput()->GetNumberOfCells() << " cells\n";

    return EXIT_SUCCESS;
}


More information about the vtk-developers mailing list