vtksmoothpolydatafilter2.cpp

   1 /*=========================================================================
   2
   3   Program:   Visualization Toolkit
   4   Module:    $RCSfile: vtksmoothpolydatafilter2.cpp,v $
   5
   6   Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
   7   All rights reserved.
   8   See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
   9
  10      This software is distributed WITHOUT ANY WARRANTY; without even
  11      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  12      PURPOSE.  See the above copyright notice for more information.
  13
  14 =========================================================================*/
  15 #include "vtksmoothpolydatafilter2.h"
  16
  17 #include "vtkCellArray.h"
  18 #include "vtkCellData.h"
  19 #include "vtkCellLocator.h"
  20 #include "vtkFloatArray.h"
  21 #include "vtkMath.h"
  22 #include "vtkInformation.h"
  23 #include "vtkInformationVector.h"
  24 #include "vtkObjectFactory.h"
  25 #include "vtkPointData.h"
  26 #include "vtkPolyData.h"
  27 #include "vtkPolygon.h"
  28 #include "vtkStreamingDemandDrivenPipeline.h"
  29 #include "vtkTriangleFilter.h"
  30 #include <QString>
  31 #include <QTextStream>
  32 #include "egvtkobject.h"
  33
  34 vtkCxxRevisionMacro(vtkSmoothPolyDataFilter2, "$Revision: 1.1.2.2 $");
  35 vtkStandardNewMacro(vtkSmoothPolyDataFilter2);
  36
  37 // The following code defines a helper class for performing mesh smoothing
  38 // across the surface of another mesh.
  39 typedef struct _vtkSmoothPoint {
  40     vtkIdType     cellId;  // cell
  41     int     subId;   // cell sub id
  42     double   p[3];    // parametric coords in cell
  43 } vtkSmoothPoint;
  44
  45 class vtkSmoothPoints { //;prevent man page generation
  46 public:
  47   vtkSmoothPoints();
  48   ~vtkSmoothPoints()
  49     {
  50     if (this->Array)
  51       {
  52       delete [] this->Array;
  53       }
  54     };
  55   vtkIdType GetNumberOfPoints() {return this->MaxId + 1;};
  56   vtkSmoothPoint *GetSmoothPoint(vtkIdType i) {return this->Array + i;};
  57   vtkSmoothPoint *InsertSmoothPoint(vtkIdType ptId)
  58     {
  59     if ( ptId >= this->Size )
  60       {
  61       this->Resize(ptId+1);
  62       }
  63     if ( ptId > this->MaxId )
  64       {
  65       this->MaxId = ptId;
  66       }
  67     return this->Array + ptId;
  68     }
  69   vtkSmoothPoint *Resize(vtkIdType sz); //reallocates data
  70   void Reset() {this->MaxId = -1;};
  71
  72   vtkSmoothPoint *Array;  // pointer to data
  73   vtkIdType MaxId;              // maximum index inserted thus far
  74   vtkIdType Size;               // allocated size of data
  75   vtkIdType Extend;             // grow array by this amount
  76 };
  77
  78 vtkSmoothPoints::vtkSmoothPoints()
  79 {
  80   this->MaxId = -1;
  81   this->Array = new vtkSmoothPoint[1000];
  82   this->Size = 1000;
  83   this->Extend = 5000;
  84 }
  85
  86 vtkSmoothPoint *vtkSmoothPoints::Resize(vtkIdType sz)
  87 {
  88   vtkSmoothPoint *newArray;
  89   vtkIdType newSize;
  90
  91   if (sz >= this->Size)
  92     {
  93     newSize = this->Size +
  94       this->Extend*(((sz-this->Size)/this->Extend)+1);
  95     }
  96   else
  97     {
  98     newSize = sz;
  99     }
 100
 101   newArray = new vtkSmoothPoint[newSize];
 102
 103   memcpy(newArray, this->Array,
 104          (sz < this->Size ? sz : this->Size) * sizeof(vtkSmoothPoint));
 105
 106   this->Size = newSize;
 107   delete [] this->Array;
 108   this->Array = newArray;
 109
 110   return this->Array;
 111 }
 112
 113 // The following code defines methods for the vtkSmoothPolyDataFilter2 class
 114 //
 115
 116 // Construct object with number of iterations 20; relaxation factor .01;
 117 // feature edge smoothing turned off; feature
 118 // angle 45 degrees; edge angle 15 degrees; and boundary smoothing turned
 119 // on. Error scalars and vectors are not generated (by default). The
 120 // convergence criterion is 0.0 of the bounding box diagonal.
 121 vtkSmoothPolyDataFilter2::vtkSmoothPolyDataFilter2()
 122 {
 123   this->Convergence = 0.0; //goes to number of specied iterations
 124   this->NumberOfIterations = 20;
 125
 126   this->RelaxationFactor = .01;
 127
 128   this->FeatureAngle = 45.0;
 129   this->EdgeAngle = 15.0;
 130   this->FeatureEdgeSmoothing = 0;
 131   this->BoundarySmoothing = 1;
 132
 133   this->GenerateErrorScalars = 0;
 134   this->GenerateErrorVectors = 0;
 135
 136   // optional second input
 137   this->SetNumberOfInputPorts(2);
 138 }
 139
 140 void vtkSmoothPolyDataFilter2::SetSource(vtkPolyData *source)
 141 {
 142   this->SetInput(1, source);
 143 }
 144
 145 vtkPolyData *vtkSmoothPolyDataFilter2::GetSource()
 146 {
 147   if (this->GetNumberOfInputConnections(1) < 1)
 148     {
 149     return NULL;
 150     }
 151   return vtkPolyData::SafeDownCast(
 152     this->GetExecutive()->GetInputData(1, 0));
 153 }
 154
 155 // #define VTK_SIMPLE_VERTEX 0
 156 // #define VTK_FIXED_VERTEX 1
 157 // #define VTK_FEATURE_EDGE_VERTEX 2
 158 // #define VTK_BOUNDARY_EDGE_VERTEX 3
 159
 160 // Special structure for marking vertices
 161 typedef struct _vtkMeshVertex
 162   {
 163   char      type;
 164   vtkIdList *edges; // connected edges (list of connected point ids)
 165 } vtkMeshVertex, *vtkMeshVertexPtr;
 166
 167 int vtkSmoothPolyDataFilter2::RequestData(
 168   vtkInformation *vtkNotUsed(request),
 169   vtkInformationVector **inputVector,
 170   vtkInformationVector *outputVector)
 171 {
 172         cout<<"=================="<<endl;
 173
 174   // get the info objects
 175   vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
 176   vtkInformation *sourceInfo = inputVector[1]->GetInformationObject(0);
 177   vtkInformation *outInfo = outputVector->GetInformationObject(0);
 178
 179   // get the input and ouptut
 180   vtkPolyData *input = vtkPolyData::SafeDownCast(
 181     inInfo->Get(vtkDataObject::DATA_OBJECT()));
 182   vtkPolyData *source = 0;
 183   if (sourceInfo)
 184     {
 185     source = vtkPolyData::SafeDownCast(
 186       sourceInfo->Get(vtkDataObject::DATA_OBJECT()));
 187     }
 188   vtkPolyData *output = vtkPolyData::SafeDownCast(
 189     outInfo->Get(vtkDataObject::DATA_OBJECT()));
 190
 191   vtkIdType numPts, numCells, i, numPolys, numStrips;
 192   int j, k;
 193   vtkIdType npts = 0;
 194   vtkIdType *pts = 0;
 195   vtkIdType p1, p2;
 196   double x[3], y[3], deltaX[3], xNew[3], conv, maxDist, dist, factor;
 197   double x1[3], x2[3], x3[3], l1[3], l2[3];
 198   double CosFeatureAngle; //Cosine of angle between adjacent polys
 199   double CosEdgeAngle; // Cosine of angle between adjacent edges
 200   double closestPt[3], dist2, *w = NULL;
 201   int iterationNumber, abortExecute;
 202   vtkIdType numSimple=0, numBEdges=0, numFixed=0, numFEdges=0;
 203   vtkPolyData *inMesh, *Mesh;
 204   vtkPoints *inPts;
 205   vtkTriangleFilter *toTris=NULL;
 206   vtkCellArray *inVerts, *inLines, *inPolys, *inStrips;
 207   vtkPoints *newPts;
 208   vtkMeshVertexPtr Verts;
 209   vtkCellLocator *cellLocator=NULL;
 210
 211   // Check input
 212   //
 213   numPts=input->GetNumberOfPoints();
 214   numCells=input->GetNumberOfCells();
 215   if (numPts < 1 || numCells < 1)
 216     {
 217     vtkErrorMacro(<<"No data to smooth!");
 218     return 1;
 219     }
 220
 221   CosFeatureAngle =
 222     cos((double) vtkMath::DegreesToRadians() * this->FeatureAngle);
 223   CosEdgeAngle = cos((double) vtkMath::DegreesToRadians() * this->EdgeAngle);
 224
 225   vtkDebugMacro(<<"Smoothing " << numPts << " vertices, " << numCells
 226                << " cells with:\n"
 227                << "\tConvergence= " << this->Convergence << "\n"
 228                << "\tIterations= " << this->NumberOfIterations << "\n"
 229                << "\tRelaxation Factor= " << this->RelaxationFactor << "\n"
 230                << "\tEdge Angle= " << this->EdgeAngle << "\n"
 231                << "\tBoundary Smoothing " << (this->BoundarySmoothing ? "On\n" : "Off\n")
 232                << "\tFeature Edge Smoothing " << (this->FeatureEdgeSmoothing ? "On\n" : "Off\n")
 233                << "\tError Scalars " << (this->GenerateErrorScalars ? "On\n" : "Off\n")
 234                << "\tError Vectors " << (this->GenerateErrorVectors ? "On\n" : "Off\n"));
 235
 236   if ( this->NumberOfIterations <= 0 || this->RelaxationFactor == 0.0)
 237     { //don't do anything! pass data through
 238     output->CopyStructure(input);
 239     output->GetPointData()->PassData(input->GetPointData());
 240     output->GetCellData()->PassData(input->GetCellData());
 241     return 1;
 242     }
 243
 244   // Peform topological analysis. What we're gonna do is build a connectivity
 245   // array of connected vertices. The outcome will be one of three
 246   // classifications for a vertex: VTK_SIMPLE_VERTEX, VTK_FIXED_VERTEX. or
 247   // VTK_EDGE_VERTEX. Simple vertices are smoothed using all connected
 248   // vertices. FIXED vertices are never smoothed. Edge vertices are smoothed
 249   // using a subset of the attached vertices.
 250   //
 251   cout<<"Analyzing topology..."<<endl;
 252         cout<<"0:numPts="<<numPts<<endl;
 253   Verts = new vtkMeshVertex[numPts];
 254   for (i=0; i<numPts; i++)
 255     {
 256         cout<<"0:VTK_SIMPLE_VERTEX"<<endl;
 257     Verts[i].type = VTK_SIMPLE_VERTEX; //can smooth
 258     Verts[i].edges = NULL;
 259     }
 260
 261   inPts = input->GetPoints();
 262   conv = this->Convergence * input->GetLength();
 263
 264   // check vertices first. Vertices are never smoothed_--------------
 265   for (inVerts=input->GetVerts(), inVerts->InitTraversal();
 266   inVerts->GetNextCell(npts,pts); )
 267     {
 268     for (j=0; j<npts; j++)
 269       {
 270         cout<<"pts[j]="<<pts[j]<<"->vertices:VTK_FIXED_VERTEX"<<endl;
 271       Verts[pts[j]].type = VTK_FIXED_VERTEX;
 272       }
 273     }
 274   this->UpdateProgress(0.10);
 275
 276   // now check lines. Only manifold lines can be smoothed------------
 277   for (inLines=input->GetLines(), inLines->InitTraversal();
 278   inLines->GetNextCell(npts,pts); )
 279     {
 280     for (j=0; j<npts; j++)
 281       {
 282         cout<<"pts[j]="<<pts[j]<<"->lines"<<endl;
 283       if ( Verts[pts[j]].type == VTK_SIMPLE_VERTEX )
 284         {
 285         if ( j == (npts-1) ) //end-of-line marked FIXED
 286           {
 287         cout<<"pts[j]="<<pts[j]<<"2:VTK_FIXED_VERTEX"<<endl;
 288           Verts[pts[j]].type = VTK_FIXED_VERTEX;
 289           }
 290         else if ( j == 0 ) //beginning-of-line marked FIXED
 291           {
 292         cout<<"pts[j]="<<pts[j]<<"3:VTK_FIXED_VERTEX"<<endl;
 293           Verts[pts[0]].type = VTK_FIXED_VERTEX;
 294           inPts->GetPoint(pts[0],x2);
 295           inPts->GetPoint(pts[1],x3);
 296           }
 297         else //is edge vertex (unless already edge vertex!)
 298           {
 299         cout<<"pts[j]="<<pts[j]<<"4:VTK_FEATURE_EDGE_VERTEX"<<endl;
 300           Verts[pts[j]].type = VTK_FEATURE_EDGE_VERTEX;
 301           Verts[pts[j]].edges = vtkIdList::New();
 302           Verts[pts[j]].edges->SetNumberOfIds(2);
 303           Verts[pts[j]].edges->SetId(0,pts[j-1]);
 304           Verts[pts[j]].edges->SetId(1,pts[j+1]);
 305           }
 306         } //if simple vertex
 307
 308       else if ( Verts[pts[j]].type == VTK_FEATURE_EDGE_VERTEX )
 309         { //multiply connected, becomes fixed!
 310         cout<<"pts[j]="<<pts[j]<<"5:VTK_FIXED_VERTEX"<<endl;
 311         Verts[pts[j]].type = VTK_FIXED_VERTEX;
 312         Verts[pts[j]].edges->Delete();
 313         Verts[pts[j]].edges = NULL;
 314         }
 315
 316       } //for all points in this line
 317     } //for all lines
 318   this->UpdateProgress(0.25);
 319
 320   // now polygons and triangle strips-------------------------------
 321   inPolys=input->GetPolys();
 322   numPolys = inPolys->GetNumberOfCells();
 323   inStrips=input->GetStrips();
 324   numStrips = inStrips->GetNumberOfCells();
 325
 326         cout<<"numPolys="<<numPolys<<endl;
 327         cout<<"numStrips="<<numStrips<<endl;
 328
 329
 330   if ( numPolys > 0 || numStrips > 0 )
 331     { //build cell structure
 332     vtkCellArray *polys;
 333     vtkIdType cellId;
 334     int numNei, nei, edge;
 335     vtkIdType numNeiPts;
 336     vtkIdType *neiPts;
 337     double normal[3], neiNormal[3];
 338     vtkIdList *neighbors;
 339
 340     neighbors = vtkIdList::New();
 341     neighbors->Allocate(VTK_CELL_SIZE);
 342
 343     inMesh = vtkPolyData::New();
 344     inMesh->SetPoints(inPts);
 345     inMesh->SetPolys(inPolys);
 346     Mesh = inMesh;
 347
 348     if ( (numStrips = inStrips->GetNumberOfCells()) > 0 )
 349       { // convert data to triangles
 350       inMesh->SetStrips(inStrips);
 351       toTris = vtkTriangleFilter::New();
 352       toTris->SetInput(inMesh);
 353       toTris->Update();
 354       Mesh = toTris->GetOutput();
 355       }
 356
 357     Mesh->BuildLinks(); //to do neighborhood searching
 358     polys = Mesh->GetPolys();
 359     this->UpdateProgress(0.375);
 360
 361     for (cellId=0, polys->InitTraversal(); polys->GetNextCell(npts,pts);
 362     cellId++)
 363       {
 364         cout<<"->cellId="<<cellId<<endl;
 365       for (i=0; i < npts; i++)
 366         {
 367           cout<<"-->i="<<i<<endl;
 368           p1 = pts[i];
 369         p2 = pts[(i+1)%npts];
 370
 371         if ( Verts[p1].edges == NULL )
 372           {
 373           Verts[p1].edges = vtkIdList::New();
 374           Verts[p1].edges->Allocate(16,6);
 375           }
 376         if ( Verts[p2].edges == NULL )
 377           {
 378           Verts[p2].edges = vtkIdList::New();
 379           Verts[p2].edges->Allocate(16,6);
 380           }
 381
 382         Mesh->GetCellEdgeNeighbors(cellId,p1,p2,neighbors);
 383         numNei = neighbors->GetNumberOfIds();
 384           cout<<"-->numNei="<<numNei<<endl;
 385
 386         edge = VTK_SIMPLE_VERTEX;
 387         if ( numNei == 0 )
 388           {
 389           edge = VTK_BOUNDARY_EDGE_VERTEX;
 390           }
 391
 392         else if ( numNei >= 2 )
 393           {
 394           // check to make sure that this edge hasn't been marked already
 395           for (j=0; j < numNei; j++)
 396             {
 397             if ( neighbors->GetId(j) < cellId )
 398               {
 399               break;
 400               }
 401             }
 402           if ( j >= numNei )
 403             {
 404             edge = VTK_FEATURE_EDGE_VERTEX;
 405             }
 406           }
 407
 408         else if ( numNei == 1 && (nei=neighbors->GetId(0)) > cellId )
 409           {
 410           vtkPolygon::ComputeNormal(inPts,npts,pts,normal);
 411           Mesh->GetCellPoints(nei,numNeiPts,neiPts);
 412           vtkPolygon::ComputeNormal(inPts,numNeiPts,neiPts,neiNormal);
 413
 414           if ( this->FeatureEdgeSmoothing &&
 415           vtkMath::Dot(normal,neiNormal) <= CosFeatureAngle )
 416             {
 417             edge = VTK_FEATURE_EDGE_VERTEX;
 418             }
 419           }
 420         else // a visited edge; skip rest of analysis
 421           {
 422           continue;
 423           }
 424
 425         if ( edge && Verts[p1].type == VTK_SIMPLE_VERTEX )
 426           {
 427           Verts[p1].edges->Reset();
 428           Verts[p1].edges->InsertNextId(p2);
 429           Verts[p1].type = edge;
 430           }
 431         else if ( (edge && Verts[p1].type == VTK_BOUNDARY_EDGE_VERTEX) ||
 432         (edge && Verts[p1].type == VTK_FEATURE_EDGE_VERTEX) ||
 433         (!edge && Verts[p1].type == VTK_SIMPLE_VERTEX ) )
 434           {
 435           Verts[p1].edges->InsertNextId(p2);
 436           if ( Verts[p1].type && edge == VTK_BOUNDARY_EDGE_VERTEX )
 437             {
 438             Verts[p1].type = VTK_BOUNDARY_EDGE_VERTEX;
 439             }
 440           }
 441
 442         if ( edge && Verts[p2].type == VTK_SIMPLE_VERTEX )
 443           {
 444           Verts[p2].edges->Reset();
 445           Verts[p2].edges->InsertNextId(p1);
 446           Verts[p2].type = edge;
 447           }
 448         else if ( (edge && Verts[p2].type == VTK_BOUNDARY_EDGE_VERTEX ) ||
 449         (edge && Verts[p2].type == VTK_FEATURE_EDGE_VERTEX) ||
 450         (!edge && Verts[p2].type == VTK_SIMPLE_VERTEX ) )
 451           {
 452           Verts[p2].edges->InsertNextId(p1);
 453           if ( Verts[p2].type && edge == VTK_BOUNDARY_EDGE_VERTEX )
 454             {
 455             Verts[p2].type = VTK_BOUNDARY_EDGE_VERTEX;
 456             }
 457           }
 458         }
 459       }
 460
 461     inMesh->Delete();
 462     if (toTris) {toTris->Delete();}
 463
 464     neighbors->Delete();
 465     }//if strips or polys
 466
 467   this->UpdateProgress(0.50);
 468
 469   //post-process edge vertices to make sure we can smooth them
 470   for (i=0; i<numPts; i++)
 471     {
 472     if ( Verts[i].type == VTK_SIMPLE_VERTEX )
 473       {
 474       numSimple++;
 475       }
 476
 477     else if ( Verts[i].type == VTK_FIXED_VERTEX )
 478       {
 479       numFixed++;
 480       }
 481
 482     else if ( Verts[i].type == VTK_FEATURE_EDGE_VERTEX ||
 483     Verts[i].type == VTK_BOUNDARY_EDGE_VERTEX )
 484       { //see how many edges; if two, what the angle is
 485
 486       if ( !this->BoundarySmoothing &&
 487       Verts[i].type == VTK_BOUNDARY_EDGE_VERTEX )
 488         {
 489         Verts[i].type = VTK_FIXED_VERTEX;
 490         numBEdges++;
 491         }
 492
 493       else if ( (npts = Verts[i].edges->GetNumberOfIds()) != 2 )
 494         {
 495         Verts[i].type = VTK_FIXED_VERTEX;
 496         numFixed++;
 497         }
 498
 499       else //check angle between edges
 500         {
 501         inPts->GetPoint(Verts[i].edges->GetId(0),x1);
 502         inPts->GetPoint(i,x2);
 503         inPts->GetPoint(Verts[i].edges->GetId(1),x3);
 504
 505         for (k=0; k<3; k++)
 506           {
 507           l1[k] = x2[k] - x1[k];
 508           l2[k] = x3[k] - x2[k];
 509           }
 510         if ( vtkMath::Normalize(l1) >= 0.0 &&
 511              vtkMath::Normalize(l2) >= 0.0 &&
 512              vtkMath::Dot(l1,l2) < CosEdgeAngle)
 513           {
 514           numFixed++;
 515           Verts[i].type = VTK_FIXED_VERTEX;
 516           }
 517         else
 518           {
 519           if ( Verts[i].type == VTK_FEATURE_EDGE_VERTEX )
 520             {
 521             numFEdges++;
 522             }
 523           else
 524             {
 525             numBEdges++;
 526             }
 527           }
 528         }//if along edge
 529       }//if edge vertex
 530     }//for all points
 531
 532   cout<<"Found\n\t" << numSimple << " simple vertices\n\t"
 533                 << numFEdges << " feature edge vertices\n\t"
 534                 << numBEdges << " boundary edge vertices\n\t"
 535                 << numFixed << " fixed vertices\n\t"<<endl;
 536         cout<<"1:numPts="<<numPts<<endl;
 537
 538     for (i=0; i<numPts; i++)
 539       {
 540         cout<<"Verts["<<i<<"].type="<<vertex_type(Verts[i].type)<<endl;
 541         if(Verts[i].edges != NULL && (npts = Verts[i].edges->GetNumberOfIds()) > 0)
 542         {
 543                 for (j=0; j<npts; j++)
 544                 {
 545                         cout<<"Verts["<<i<<"].edges->GetId("<<j<<")="<<Verts[i].edges->GetId(j)<<endl;
 546                 }//for all connected points
 547         }
 548       }
 549
 550   cout<<"Beginning smoothing iterations..."<<endl;
 551         cout<<"2:numPts="<<numPts<<endl;
 552
 553   // We've setup the topology...now perform Laplacian smoothing
 554   //
 555   newPts = vtkPoints::New();
 556   newPts->SetNumberOfPoints(numPts);
 557
 558   // If Source defined, we do constrained smoothing (that is, points are
 559   // constrained to the surface of the mesh object).
 560   if ( source )
 561     {
 562     this->SmoothPoints = new vtkSmoothPoints;
 563     vtkSmoothPoint *sPtr;
 564     cellLocator = vtkCellLocator::New();
 565     w = new double[source->GetMaxCellSize()];
 566
 567     cellLocator->SetDataSet(source);
 568     cellLocator->BuildLocator();
 569
 570     for (i=0; i < numPts; i++)
 571       {
 572       sPtr = this->SmoothPoints->InsertSmoothPoint(i);
 573       cellLocator->FindClosestPoint(inPts->GetPoint(i), closestPt,
 574                                     sPtr->cellId, sPtr->subId, dist2);
 575       newPts->SetPoint(i, closestPt);
 576       }
 577     }
 578   else //smooth normally
 579     {
 580     for (i=0; i<numPts; i++) //initialize to old coordinates
 581       {
 582       newPts->SetPoint(i,inPts->GetPoint(i));
 583       }
 584     }
 585
 586   factor = this->RelaxationFactor;
 587   for ( maxDist=VTK_DOUBLE_MAX, iterationNumber=0, abortExecute=0;
 588   maxDist > conv && iterationNumber < this->NumberOfIterations && !abortExecute;
 589   iterationNumber++ )
 590     {
 591
 592     if ( iterationNumber && !(iterationNumber % 5) )
 593       {
 594       this->UpdateProgress (0.5 + 0.5*iterationNumber/this->NumberOfIterations);
 595       if (this->GetAbortExecute())
 596         {
 597         cout<<"ABOOOOOOOOOOOOOOOORT!!!!!!!!!!"<<endl;
 598         abortExecute = 1;
 599         break;
 600         }
 601       }
 602
 603     maxDist=0.0;
 604 //      cout<<"numPts="<<numPts<<endl;
 605     for (i=0; i<numPts; i++)
 606       {
 607 //              cout<<"Verts["<<i<<"].type="<<(int)Verts[i].type<<endl;
 608       if ( Verts[i].type != VTK_FIXED_VERTEX && Verts[i].edges != NULL &&
 609       (npts = Verts[i].edges->GetNumberOfIds()) > 0 )
 610         {
 611         newPts->GetPoint(i, x); //use current points
 612         deltaX[0] = deltaX[1] = deltaX[2] = 0.0;
 613         for (j=0; j<npts; j++)
 614           {
 615           newPts->GetPoint(Verts[i].edges->GetId(j), y);
 616           for (k=0; k<3; k++)
 617             {
 618             deltaX[k] += (y[k] - x[k]) / npts;
 619             }
 620           }//for all connected points
 621
 622         for (k=0;k<3;k++)
 623           {
 624           xNew[k] = x[k] + factor * deltaX[k];
 625           }
 626
 627         // Constrain point to surface
 628         if ( source )
 629           {
 630           vtkSmoothPoint *sPtr = this->SmoothPoints->GetSmoothPoint(i);
 631           vtkCell *cell=NULL;
 632
 633           if ( sPtr->cellId >= 0 ) //in cell
 634             {
 635             cell = source->GetCell(sPtr->cellId);
 636             }
 637
 638           if ( !cell || cell->EvaluatePosition(xNew, closestPt,
 639           sPtr->subId, sPtr->p, dist2, w) == 0)
 640             { // not in cell anymore
 641             cellLocator->FindClosestPoint(xNew, closestPt, sPtr->cellId,
 642                                           sPtr->subId, dist2);
 643             }
 644           for (k=0; k<3; k++)
 645             {
 646             xNew[k] = closestPt[k];
 647             }
 648           }
 649
 650         newPts->SetPoint(i,xNew);
 651         if ( (dist = vtkMath::Norm(deltaX)) > maxDist )
 652           {
 653           maxDist = dist;
 654           }
 655         }//if can move point
 656       }//for all points
 657     } //for not converged or within iteration count
 658
 659   vtkDebugMacro(<<"Performed " << iterationNumber << " smoothing passes");
 660   if ( source )
 661     {
 662     cellLocator->Delete();
 663     delete this->SmoothPoints;
 664     delete [] w;
 665     }
 666
 667   // Update output. Only point coordinates have changed.
 668   //
 669   output->GetPointData()->PassData(input->GetPointData());
 670   output->GetCellData()->PassData(input->GetCellData());
 671
 672   if ( this->GenerateErrorScalars )
 673     {
 674     vtkFloatArray *newScalars = vtkFloatArray::New();
 675     newScalars->SetNumberOfTuples(numPts);
 676     for (i=0; i<numPts; i++)
 677       {
 678       inPts->GetPoint(i,x1);
 679       newPts->GetPoint(i,x2);
 680       newScalars->SetComponent(i,0,
 681                                sqrt(vtkMath::Distance2BetweenPoints(x1,x2)));
 682         cout<<"dist["<<i<<"]="<<sqrt(vtkMath::Distance2BetweenPoints(x1,x2))<<endl;
 683       }
 684     int idx = output->GetPointData()->AddArray(newScalars);
 685         cout<<"ERROR SCALARS STORED IN ARRAY "<<idx<<endl;
 686     output->GetPointData()->SetActiveAttribute(idx, vtkDataSetAttributes::SCALARS);
 687     for (i=0; i<numPts; i++)
 688     {
 689       double dist=newScalars->GetComponent(i-1,idx);
 690       cout<<"dist["<<i<<"]="<<dist<<endl;
 691     }
 692
 693     newScalars->Delete();
 694     }
 695
 696   if ( this->GenerateErrorVectors )
 697     {
 698     vtkFloatArray *newVectors = vtkFloatArray::New();
 699     newVectors->SetNumberOfComponents(3);
 700     newVectors->SetNumberOfTuples(numPts);
 701     for (i=0; i<numPts; i++)
 702       {
 703       inPts->GetPoint(i,x1);
 704       newPts->GetPoint(i,x2);
 705       for (j=0; j<3; j++)
 706         {
 707         x3[j] = x2[j] - x1[j];
 708         }
 709       newVectors->SetTuple(i,x3);
 710       }
 711     output->GetPointData()->SetVectors(newVectors);
 712     newVectors->Delete();
 713     }
 714
 715   output->SetPoints(newPts);
 716   newPts->Delete();
 717
 718   output->SetVerts(input->GetVerts());
 719   output->SetLines(input->GetLines());
 720   output->SetPolys(input->GetPolys());
 721   output->SetStrips(input->GetStrips());
 722
 723   //free up connectivity storage
 724   for (i=0; i<numPts; i++)
 725     {
 726     if ( Verts[i].edges != NULL )
 727       {
 728       Verts[i].edges->Delete();
 729       Verts[i].edges = NULL;
 730       }
 731     }
 732   delete [] Verts;
 733
 734   return 1;
 735 }
 736
 737 int vtkSmoothPolyDataFilter2::FillInputPortInformation(int port,
 738                                                       vtkInformation *info)
 739 {
 740   if (!this->Superclass::FillInputPortInformation(port, info))
 741     {
 742     return 0;
 743     }
 744
 745   if (port == 1)
 746     {
 747     info->Set(vtkAlgorithm::INPUT_IS_OPTIONAL(), 1);
 748     }
 749   return 1;
 750 }
 751
 752 void vtkSmoothPolyDataFilter2::PrintSelf(ostream& os, vtkIndent indent)
 753 {
 754   this->Superclass::PrintSelf(os,indent);
 755
 756   os << indent << "Convergence: " << this->Convergence << "\n";
 757   os << indent << "Number of Iterations: " << this->NumberOfIterations << "\n";
 758   os << indent << "Relaxation Factor: " << this->RelaxationFactor << "\n";
 759   os << indent << "Feature Edge Smoothing: " << (this->FeatureEdgeSmoothing ? "On\n" : "Off\n");
 760   os << indent << "Feature Angle: " << this->FeatureAngle << "\n";
 761   os << indent << "Edge Angle: " << this->EdgeAngle << "\n";
 762   os << indent << "Boundary Smoothing: " << (this->BoundarySmoothing ? "On\n" : "Off\n");
 763   os << indent << "Generate Error Scalars: " << (this->GenerateErrorScalars ? "On\n" : "Off\n");
 764   os << indent << "Generate Error Vectors: " << (this->GenerateErrorVectors ? "On\n" : "Off\n");
 765   if ( this->GetSource() )
 766     {
 767     os << indent << "Source: " << (void *)this->GetSource() << "\n";
 768     }
 769   else
 770     {
 771     os << indent << "Source (none)\n";
 772     }
 773
 774 }