vtkConnectedTubeFilter.C 25.3 KB
Newer Older
hrchilds's avatar
hrchilds committed
1 2
/*****************************************************************************
*
bonnell's avatar
bonnell committed
3
* Copyright (c) 2000 - 2017, Lawrence Livermore National Security, LLC
hrchilds's avatar
hrchilds committed
4
* Produced at the Lawrence Livermore National Laboratory
5
* LLNL-CODE-442911
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 "vtkConnectedTubeFilter.h"
hrchilds's avatar
hrchilds committed
40 41 42 43 44

#include <vtkCellArray.h>
#include <vtkCellData.h>
#include <vtkCleanPolyData.h>
#include <vtkFloatArray.h>
45 46
#include <vtkInformation.h>
#include <vtkInformationVector.h>
hrchilds's avatar
hrchilds committed
47 48 49
#include <vtkMath.h>
#include <vtkObjectFactory.h>
#include <vtkPointData.h>
hrchilds's avatar
hrchilds committed
50
#include <vtkPolyData.h>
hrchilds's avatar
hrchilds committed
51 52
#include <vtkPolyLine.h>
#include <vtkTubeFilter.h>
hrchilds's avatar
hrchilds committed
53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 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
// ----------------------------------------------------------------------------
//                            class PointSequence
// ----------------------------------------------------------------------------


// ****************************************************************************
//  Constructor:  PointSequence::PointSequence
//
//  Programmer:  Jeremy Meredith
//  Creation:    November  1, 2002
//
// ****************************************************************************
vtkConnectedTubeFilter::PointSequence::PointSequence()
{
    length    = 0;
    index     = NULL;
    cellindex = NULL;
}

// ****************************************************************************
//  Denstructor:  PointSequence::~PointSequence
//
//  Programmer:  Jeremy Meredith
//  Creation:    November  1, 2002
//
// ****************************************************************************
vtkConnectedTubeFilter::PointSequence::~PointSequence()
{
    if (index)
        delete[] index;
    index = NULL;
    if (cellindex)
        delete[] cellindex;
    cellindex = NULL;
}

// ****************************************************************************
//  Method:  PointSequence::Init
//
//  Purpose:
//    Initialize the data structures now that we know how big it might be
//
//  Programmer:  Jeremy Meredith
//  Creation:    November  1, 2002
//
//  Modifications:
//    Hank Childs (for Jeremy Meredith), Mon Apr  7 11:46:51 PDT 2003
//    Increased size of cellindex.
//
102 103 104
//
//  Jean Favre, Tue May  7 16:38:37 CEST 2013
//  Used vtkIdType where needed
hrchilds's avatar
hrchilds committed
105 106 107 108 109 110 111 112 113 114
// ****************************************************************************
void
vtkConnectedTubeFilter::PointSequence::Init(int maxlen)
{
    if (index)
        delete[] index;
    if (cellindex)
        delete[] cellindex;

    length = 0;
115 116
    index  = new vtkIdType[maxlen];
    cellindex = new vtkIdType[maxlen];
hrchilds's avatar
hrchilds committed
117 118 119 120 121 122 123 124 125 126 127
}

// ****************************************************************************
//  Method:  PointSequence::Add
//
//  Purpose:
//    Adds a new point, with its cell index, to the current sequence
//
//  Programmer:  Jeremy Meredith
//  Creation:    November  1, 2002
//
128 129
//  Jean Favre, Tue May  7 16:38:37 CEST 2013
//  Used vtkIdType where needed
hrchilds's avatar
hrchilds committed
130 131
// ****************************************************************************
void
132
vtkConnectedTubeFilter::PointSequence::Add(vtkIdType i, vtkIdType ci)
hrchilds's avatar
hrchilds committed
133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
{
    index[length]     = i;
    cellindex[length] = ci;
    length++;
}


// ----------------------------------------------------------------------------
//                          class PointSequenceList
// ----------------------------------------------------------------------------


// ****************************************************************************
//  Constructor:  PointSequenceList::PointSequenceList
//
//  Programmer:  Jeremy Meredith
//  Creation:    November  1, 2002
//
151 152 153 154
//  Modifications:
//    Rich Cook and Hank Childs, Thu Oct  2 16:32:55 PDT 2008
//    Initialized data member used for supporting loops.
//
hrchilds's avatar
hrchilds committed
155 156 157 158 159 160 161 162 163 164 165 166
// ****************************************************************************
vtkConnectedTubeFilter::PointSequenceList::PointSequenceList()
{
    pts             = NULL;
    len             = 0;
    numneighbors    = NULL;
    connectivity[0] = NULL;
    connectivity[1] = NULL;
    cellindex       = NULL;

    visited = NULL;
    index   = -1;
167
    lookforloops = false;
hrchilds's avatar
hrchilds committed
168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
}

// ****************************************************************************
//  Denstructor:  PointSequenceList::~PointSequenceList
//
//  Programmer:  Jeremy Meredith
//  Creation:    November  1, 2002
//
// ****************************************************************************
vtkConnectedTubeFilter::PointSequenceList::~PointSequenceList()
{
    if (visited)
        delete[] visited;
    if (numneighbors)
        delete[] numneighbors;
    if (connectivity[0])
        delete[] connectivity[0];
    if (connectivity[1])
        delete[] connectivity[1];
    if (cellindex)
        delete[] cellindex;
    // Don't delete 'pts', because we don't own it.
}

// ****************************************************************************
//  Method:  PointSequenceList::Build
//
//  Purpose:
//    Extract the 
//
//  Programmer:  Jeremy Meredith
//  Creation:    November  1, 2002
//
201 202
//  Jean Favre, Tue May  7 16:38:37 CEST 2013
//  Used vtkIdType where needed
hrchilds's avatar
hrchilds committed
203 204 205 206 207
// ****************************************************************************
bool 
vtkConnectedTubeFilter::PointSequenceList::Build(vtkPoints *points,
                                                 vtkCellArray *lines)
{
bonnell's avatar
bonnell committed
208
    pts             = points;
hrchilds's avatar
hrchilds committed
209
    len             = points->GetNumberOfPoints();
210 211 212 213
    numneighbors    = new vtkIdType[len];
    connectivity[0] = new vtkIdType[len];
    connectivity[1] = new vtkIdType[len];
    cellindex       = new vtkIdType[len];
hrchilds's avatar
hrchilds committed
214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233

    vtkIdType *cells = lines->GetPointer();

    // Initalize all points to be disconnected from each other
    int i;
    for (i=0; i<len; i++)
    {
        numneighbors[i] = 0;
    }

    int numCells = lines->GetNumberOfCells();
    for (i=0; i<numCells; i++)
    {
        // We assume all cells are two-point lines (i.e. not polylines)
        if (cells[i*3] != 2)
        {
            return false;
        }

        // Get the begin and end index for this segment
234 235
        vtkIdType a = cells[i*3 + 1];
        vtkIdType b = cells[i*3 + 2];
hrchilds's avatar
hrchilds committed
236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262

        // If we have two neighbors already, this is a T intersection
        if (numneighbors[a] >= 2 || numneighbors[b] >= 2)
        {
            return false;
        }

        // Set the neighbors of each endpoint to be each other
        connectivity[numneighbors[a]][a] = b;
        connectivity[numneighbors[b]][b] = a;
        numneighbors[a]++;
        numneighbors[b]++;
        cellindex[a] = i;
        cellindex[b] = i;
    }
    return true;
}

// ****************************************************************************
//  Method:  PointSequenceList::InitTraversal
//
//  Purpose:
//    Create a new 'visited' array and reset the iterator index
//
//  Programmer:  Jeremy Meredith
//  Creation:    November  1, 2002
//
263 264 265 266
//  Modifications:
//    Rich Cook and Hank Childs, Thu Oct  2 16:32:55 PDT 2008
//    Initialized data member used for supporting loops.
//
hrchilds's avatar
hrchilds committed
267 268 269 270 271 272 273 274 275 276 277 278
// ****************************************************************************
void
vtkConnectedTubeFilter::PointSequenceList::InitTraversal()
{
    if (visited)
        delete[] visited;

    visited = new bool[len];
    for (int i=0; i<len; i++)
        visited[i] = false;

    index = 0;
279
    lookforloops = false;
hrchilds's avatar
hrchilds committed
280 281 282 283 284 285 286 287 288 289 290 291 292 293
}

// ****************************************************************************
//  Method:  PointSequenceList::GetNextSequence
//
//  Purpose:
//    Extract another sequence from the connectivity array and return it
//    as a PointSequence.  Skip repeated adjacent points.
//
//  Returns:  true if we found another sequence; false otherwise
//
//  Programmer:  Jeremy Meredith
//  Creation:    November  1, 2002
//
294 295 296 297
//  Modifications:
//    Rich Cook and Hank Childs, Thu Oct  2 16:32:55 PDT 2008
//    Added support for loops.
//
298 299 300 301 302
//    Eric Brugger, Tue Dec 29 16:35:34 PST 2009
//    I modified the logic that looks for loops to only consider points
//    that have 2 neighbors.  Previously it would have considered points
//    with no neighbors, which caused it to reference uninitialized memory.
//
303 304 305 306 307
//    Jean Favre, Tue May  7 16:38:37 CEST 2013
//    I modified the calls to GetPoint() to use the other variant of the call
//    with two arguments. The previous version would never succeed in the test
//    to remove sequential identical points.
//    Used vtkIdType where needed
hrchilds's avatar
hrchilds committed
308 309 310 311 312 313 314 315
// ****************************************************************************
bool
vtkConnectedTubeFilter::PointSequenceList::GetNextSequence(PointSequence &seq)
{
    // index is already set by InitTraversal and previous calls to this fn
    for (; index < len; index++)
    {
        // if numneighbors is 1, then this is a start point
316 317
        if (((lookforloops && numneighbors[index] == 2) ||
              numneighbors[index] == 1) && !visited[index])
hrchilds's avatar
hrchilds committed
318
        {
319 320
            vtkIdType current = index;
            vtkIdType previous = -1;
hrchilds's avatar
hrchilds committed
321 322 323 324 325
            seq.Init(len);
            seq.Add(current, cellindex[current]);
            visited[current] = true;
            while (true)
            {
326 327 328
                vtkIdType n1       = connectivity[0][current];
                vtkIdType n2       = connectivity[1][current];
                vtkIdType next     = (n1 == previous) ? n2 : n1;
hrchilds's avatar
hrchilds committed
329 330 331
                previous = current;
                current = next;

332
 
hrchilds's avatar
hrchilds committed
333 334
                // we must skip any sequential identical points:
                // 1) they are useless, and 2) they mess up calculations
335 336 337 338
                double prePt[3];
                double curPt[3];
                pts->GetPoint(previous, prePt);
                pts->GetPoint(current, curPt);
bonnell's avatar
bonnell committed
339 340 341
                if (prePt[0] != curPt[0] ||
                    prePt[1] != curPt[1] ||
                    prePt[2] != curPt[2])
hrchilds's avatar
hrchilds committed
342 343 344 345
                {
                    seq.Add(current, cellindex[current]);
                }

346 347 348 349 350 351 352
                // check for a loop (AFTER adding the node again...)
                if (lookforloops && visited[current]) 
                {
                    break; 
                }
                visited[current] = true;

hrchilds's avatar
hrchilds committed
353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371
                bool endpoint = (numneighbors[current] == 1);
                if (endpoint)
                    break;
            }

            // We may have one duplicated cell index due to the way the
            // connectivity was built -- shift them all down by one from
            // the right spot, so the only duplicated one is the last one.
            for (int i=1; i<seq.length-1; i++)
            {
                if (seq.cellindex[i-1] == seq.cellindex[i])
                    seq.cellindex[i] = seq.cellindex[i+1];
            }

            // true ==> success; got another sequence
            return true;
        }
    }

372 373 374 375 376 377 378
    if (index == len && !lookforloops) 
    {
        lookforloops = true; 
        index=0; 
        return GetNextSequence(seq);
    }

hrchilds's avatar
hrchilds committed
379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401
    // false ==> failed; no more sequences
    return false;
}


// ----------------------------------------------------------------------------
//                        class vtkConnectedTubeFilter
// ----------------------------------------------------------------------------


vtkStandardNewMacro(vtkConnectedTubeFilter);

// ****************************************************************************
//  Constructor:  vtkConnectedTubeFilter::vtkConnectedTubeFilter
//
//  Programmer:  Jeremy Meredith
//  Creation:    November  1, 2002
//
// ****************************************************************************
vtkConnectedTubeFilter::vtkConnectedTubeFilter()
{
    this->Radius = 0.5;
    this->NumberOfSides = 3;
bonnell's avatar
bonnell committed
402 403
    this->CreateNormals = false;
    this->Capping = false;
hrchilds's avatar
hrchilds committed
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

    pseqlist = NULL;
}

// ****************************************************************************
//  Denstructor:  vtkConnectedTubeFilter::~vtkConnectedTubeFilter
//
//  Programmer:  Jeremy Meredith
//  Creation:    November  1, 2002
//
// ****************************************************************************
vtkConnectedTubeFilter::~vtkConnectedTubeFilter()
{
    if (pseqlist)
        delete pseqlist;
    pseqlist = NULL;
}

// ****************************************************************************
//  Method:  vtkConnectedTubeFilter::BuildConnectivityArrays
//
//  Purpose:
//    Create our new style connectivity arrays.
//
//  Programmer:  Jeremy Meredith
//  Creation:    November  1, 2002
//
431 432 433 434
//  Modifications:
//    Eric Brugger, Wed Jan  9 11:31:17 PST 2013
//    Modified to inherit from vtkPolyDataAlgorithm.
//
hrchilds's avatar
hrchilds committed
435
// ****************************************************************************
436
bool vtkConnectedTubeFilter::BuildConnectivityArrays(vtkPolyData *input)
hrchilds's avatar
hrchilds committed
437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467
{
    vtkPoints    *inPts   = NULL;
    vtkCellArray *inLines = NULL;
    int numPts;
    int numCells;
    vtkDebugMacro(<<"Building tube connectivity arrays");

    if (!(inPts=input->GetPoints())               || 
        (numPts = inPts->GetNumberOfPoints()) < 1 ||
        !(inLines = input->GetLines())            ||
        (numCells = inLines->GetNumberOfCells()) < 1)
    {
        vtkDebugMacro(<< ": No input data!\n");
        return false;
    }

    // Build the new connectivity, and check for errors while we're at it
    pseqlist = new PointSequenceList;
    bool success = pseqlist->Build(inPts, inLines);

    if (!success)
    {
        // We know we can't use it anyway; delete it now
        delete pseqlist;
        pseqlist = NULL;
    }

    return success;
}

// ****************************************************************************
468
//  Method:  vtkConnectedTubeFilter::RequestData
hrchilds's avatar
hrchilds committed
469 470 471 472 473 474 475 476 477 478 479
//
//  Purpose:
//    Normal vtk filter execution.
//
//  Programmer:  Jeremy Meredith
//  Creation:    November  1, 2002
//
//  Modifications:
//    Hank Childs, Mon Apr  7 10:02:31 PDT 2003
//    Allocate memory for points, because VTK does not do that for you.
//
480 481 482 483 484
//    Jeremy Meredith, Wed Oct 15 15:56:34 EDT 2008
//    It appears we were creating the tube points incorrectly (creating
//    degenerate polys between 360 and 0 degrees) by partitioning into
//    npts-1 sides, instead of npts sides.  Fixed that.
//
bonnell's avatar
bonnell committed
485 486 487
//    Kathleen Biagas, Thu Sep 6 11:15:29 MST 2012
//    Preserve coordinate data type.
//
488 489 490
//    Eric Brugger, Wed Jan  9 11:31:17 PST 2013
//    Modified to inherit from vtkPolyDataAlgorithm.
//
491 492 493 494
//    Jean Favre, Tue May  7 16:38:37 CEST 2013
//    I modified the calls to GetPoint() to use the other variant of the call
//    with two arguments. The previous version would fail getting the right values.
//    Used vtkIdType where needed
495 496 497 498 499
//
//    Jean Favre, Tue May 28 13:28:41 CEST 2013
//    I added a quaternion interpolation to smoothly interpolate the direction vectors
//    used to construct the tube. This eliminates abrupt twists along the tube
//
hrchilds's avatar
hrchilds committed
500
// ****************************************************************************
501 502 503 504
int vtkConnectedTubeFilter::RequestData(
    vtkInformation *vtkNotUsed(request),
    vtkInformationVector **inputVector,
    vtkInformationVector *outputVector)
hrchilds's avatar
hrchilds committed
505
{
506 507 508 509 510 511 512 513 514 515 516 517
    // get the info objects
    vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
    vtkInformation *outInfo = outputVector->GetInformationObject(0);

    //
    // Initialize some frequently used values.
    //
    vtkPolyData  *input = vtkPolyData::SafeDownCast(
        inInfo->Get(vtkDataObject::DATA_OBJECT()));
    vtkPolyData *output = vtkPolyData::SafeDownCast(
        outInfo->Get(vtkDataObject::DATA_OBJECT()));

hrchilds's avatar
hrchilds committed
518 519 520 521 522 523 524 525 526 527 528 529 530 531
    // Get all the appropriate input arrays
    vtkPoints    *inPts   = NULL;
    vtkCellArray *inLines = NULL;
    vtkCellData  *inCD    = input->GetCellData();
    vtkPointData *inPD    = input->GetPointData();
    int numPts;
    int numCells;
    vtkDebugMacro(<<"Creating tube");

    // We assume BuildConnectivityArrays was already called
    if (!pseqlist)
    {
        vtkErrorMacro(<< ": Connectivity was not built yet; need to call "
                         "vtkConnectedTubeFilter::BuildConnectivityArrays()\n");
532
        return 1;
hrchilds's avatar
hrchilds committed
533 534 535 536 537 538 539 540
    }

    if (!(inPts=input->GetPoints())               || 
        (numPts = inPts->GetNumberOfPoints()) < 1 ||
        !(inLines = input->GetLines())            ||
        (numCells = inLines->GetNumberOfCells()) < 1)
    {
        vtkDebugMacro(<< ": No input data!\n");
541
        return 1;
hrchilds's avatar
hrchilds committed
542 543 544 545 546
    }

    // Set up the output arrays
    int maxNewCells  = numCells * (NumberOfSides + 2);
    int maxNewPoints = numCells * NumberOfSides * 2;
bonnell's avatar
bonnell committed
547
    vtkPoints     *newPts     = vtkPoints::New(inPts->GetDataType());
hrchilds's avatar
hrchilds committed
548 549 550
    newPts->Allocate(maxNewPoints);
    vtkCellArray  *newCells   = vtkCellArray::New();
    newCells->Allocate(maxNewCells, 4*maxNewCells + 2*NumberOfSides);
bonnell's avatar
bonnell committed
551
    vtkDataArray *newNormals = NULL;
hrchilds's avatar
hrchilds committed
552 553 554 555 556 557 558 559 560 561 562 563

    vtkPointData  *newPD      = NULL;
    newPD = output->GetPointData();
    newPD->CopyNormalsOff();
    newPD->CopyAllocate(inPD, maxNewPoints);

    vtkCellData   *newCD      = NULL;
    newCD = output->GetCellData();
    newCD->CopyAllocate(inCD, maxNewCells);
    
    if (CreateNormals)
    {
bonnell's avatar
bonnell committed
564
        newNormals = inPts->GetData()->NewInstance();
hrchilds's avatar
hrchilds committed
565 566 567 568 569 570 571 572 573 574 575 576 577 578 579
        newNormals->SetNumberOfComponents(3);
        newNormals->SetName("Normals");
        newPD->SetNormals(newNormals);
        newNormals->Delete();
    }

    // Iterate over all connected point sequences
    PointSequence seq;
    pseqlist->InitTraversal();
    while (pseqlist->GetNextSequence(seq))
    {
        // Skip any standalone points.
        if (seq.length == 1)
            continue;

580
        double prev_dir[3] = {0., 0., 0.}, prev_v1[3], prev_v2[3];
hrchilds's avatar
hrchilds committed
581 582 583 584 585 586
        for (int i=0; i<seq.length; i++)
        {
            bool firstPoint = (i==0);
            bool lastPoint  = (i==seq.length-1);

            // Get the current, previous, and next indices
587 588 589
            vtkIdType ix  = seq.index[i];
            vtkIdType ix1 = (firstPoint ?  seq.index[i] : seq.index[i-1]);
            vtkIdType ix2 = (lastPoint  ?  seq.index[i] : seq.index[i+1]);
hrchilds's avatar
hrchilds committed
590 591

            // Use a centered difference approximation for direction
592 593 594 595
            double v0[3], v1[3], v2[3], dir[3], pt[3];
            inPts->GetPoint(ix2, v2);
            inPts->GetPoint(ix1, v1);
            vtkMath::Subtract(v2, v1, dir);
hrchilds's avatar
hrchilds committed
596 597 598 599 600 601

            // If our centered difference was zero, do a forward
            // difference instead.  We ensured no sequential points
            // are identical, so this can't fail.
            if (dir[0]==0 && dir[1]==0 && dir[2]==0)
            {
602 603
               inPts->GetPoint(ix, v0);
               vtkMath::Subtract(v2, v0, dir);
hrchilds's avatar
hrchilds committed
604
            }
605 606
            //inNormals->GetTuple(i, v1);
            vtkMath::Cross(dir, v1, v2);
hrchilds's avatar
hrchilds committed
607 608
            vtkMath::Normalize(v2);

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
            if(firstPoint)
            {
              // Get a couple vectors orthogonal to our direction
              //double v1[3], v2[3];
              vtkMath::Perpendiculars(dir, v1,v2, 0.0);
              vtkMath::Normalize(v1);
              vtkMath::Normalize(v2);
            }
            else
            {
              // we can't just get two new orthogonal vectors each time because of abrupt changes
              // We will construct the matrix which xforms the "previous dir" vector
              // into this new dir vector, and then use this rigid body motion to
              // xform v1 and v2 from the previous point
              double Xform[3][3];
              double vec[3];
              vtkMath::Cross(prev_dir, dir, vec);
              double costheta = vtkMath::Dot(prev_dir, dir);
              double sintheta = vtkMath::Norm(vec);
              double theta = atan2(sintheta, costheta);
              if (sintheta != 0)
                {
                vec[0] /= sintheta;
                vec[1] /= sintheta;
                vec[2] /= sintheta;
                }
              // convert to quaternion
              costheta = cos(0.5*theta);
              sintheta = sin(0.5*theta);
              double quat[4];
              quat[0] = costheta;
              quat[1] = vec[0]*sintheta;
              quat[2] = vec[1]*sintheta;
              quat[3] = vec[2]*sintheta;
              // convert to matrix
              vtkMath::QuaternionToMatrix3x3(quat, Xform); // Xform will xform prev_dir into dir
              // now use it
              vtkMath::Multiply3x3 (Xform, prev_v1, v1);
              vtkMath::Multiply3x3 (Xform, prev_v2, v2);
              vtkMath::Normalize(v1);
              vtkMath::Normalize(v2);
            }

            prev_dir[0] = dir[0]; prev_dir[1] = dir[1]; prev_dir[2] = dir[2];
            prev_v1[0]  = v1[0];  prev_v1[1]  = v1[1];  prev_v1[2]  = v1[2];
            prev_v2[0]  = v2[0];  prev_v2[1]  = v2[1];  prev_v2[2]  = v2[2];

hrchilds's avatar
hrchilds committed
656 657 658
            // Hang on to the first point index we create; we need it 
            // to create the cells
            vtkIdType firstIndex = newPts->GetNumberOfPoints();
bonnell's avatar
bonnell committed
659
    
660 661
            //double *pt = inPts->GetPoint(ix);
            inPts->GetPoint(ix, pt);
hrchilds's avatar
hrchilds committed
662 663
            for (int j = 0 ; j < NumberOfSides ; j++)
            {
bonnell's avatar
bonnell committed
664 665 666 667
                double q = (j * 2. * vtkMath::Pi()) / double(NumberOfSides);
                double sq = sin(q);
                double cq = cos(q);
                double normal[3] = { v1[0]*cq + v2[0]*sq,
hrchilds's avatar
hrchilds committed
668 669
                                    v1[1]*cq + v2[1]*sq,
                                    v1[2]*cq + v2[2]*sq};
bonnell's avatar
bonnell committed
670 671 672
                double x = pt[0] + Radius * normal[0];
                double y = pt[1] + Radius * normal[1];
                double z = pt[2] + Radius * normal[2];
hrchilds's avatar
hrchilds committed
673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692
                vtkIdType id = newPts->InsertNextPoint(x,y,z);
                if (CreateNormals)
                    newNormals->InsertNextTuple(normal);
                newPD->CopyData(inPD, ix, id);
            }

            // Do the start cap
            if (firstPoint && Capping)
            {
                vtkIdType id = newCells->InsertNextCell(NumberOfSides);
                for (int j=0; j<NumberOfSides; j++)
                    newCells->InsertCellPoint(firstIndex + j);
                newCD->CopyData(inCD, seq.cellindex[i], id);
            }

            // Do the next segment
            if (!firstPoint)
            {
                for (int j=0; j<NumberOfSides; j++)
                {
693
                    vtkIdType p[4] =
hrchilds's avatar
hrchilds committed
694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728
                    {
                        firstIndex + j,
                        firstIndex + j - NumberOfSides,
                        firstIndex + ((j+1) % NumberOfSides) - NumberOfSides,
                        firstIndex + ((j+1) % NumberOfSides)
                    };
                    vtkIdType id = newCells->InsertNextCell(4, p);
                    newCD->CopyData(inCD, seq.cellindex[i-1], id);
                }
            }

            // Do the end cap
            if (lastPoint && Capping)
            {
                vtkIdType id = newCells->InsertNextCell(NumberOfSides);
                for (int j=0; j<NumberOfSides; j++)
                    newCells->InsertCellPoint((firstIndex+NumberOfSides-1) - j);
                newCD->CopyData(inCD, seq.cellindex[i-1], id);
            }
        }
    }

    // Final cleanup
    newPD->Squeeze();
    newCD->Squeeze();

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

    output->SetPolys(newCells);
    newCells->Delete();

    // don't forget the sequence list; we're done with it
    delete pseqlist;
    pseqlist = NULL;
729 730

    return 1;
hrchilds's avatar
hrchilds committed
731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752
}

// ****************************************************************************
//  Method:  vtkConnectedTubeFilter::PrintSelf
//
//  Purpose:
//    Print myself.
//
//  Programmer:  Jeremy Meredith
//  Creation:    November  1, 2002
//
// ****************************************************************************
void vtkConnectedTubeFilter::PrintSelf(ostream& os, vtkIndent indent)
{
    this->Superclass::PrintSelf(os,indent);

    os << indent << "Radius: " << this->Radius << "\n";
    os << indent << "Number Of Sides: " << this->NumberOfSides << "\n";
    os << indent << "Create Normals: " 
       << (this->CreateNormals ? "On\n" : "Off\n");
    os << indent << "Capping: " << this->Capping << endl;
}