vtkSQ3SubdivisionFilter.cxx Source File

00001 /*=========================================================================
00002 
00003   Program:   Visualization Toolkit
00004   Module:    $RCSfile: vtkSQ3SubdivisionFilter.cxx,v $
00005   Language:  C++
00006   Date:      $Date: 2002/12/17 02:05:38 $
00007   Version:   $Revision: 1.20 $
00008 
00009   Copyright (c) 1993-2002 Ken Martin, Will Schroeder, Bill Lorensen 
00010   All rights reserved.
00011   See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
00012 
00013      This software is distributed WITHOUT ANY WARRANTY; without even 
00014      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
00015      PURPOSE.  See the above copyright notice for more information.
00016 
00017 =========================================================================*/
00018 #include "StdAfx.h"
00019 
00020 #include "vtkSQ3SubdivisionFilter.h"
00021 
00022 #include "vtkObjectFactory.h"
00023 #include "vtkCellArray.h"
00024 #include "vtkCellData.h"
00025 #include "vtkEdgeTable.h"
00026 #include "vtkPointData.h"
00027 #include "vtkPolyData.h"
00028 
00029 vtkCxxRevisionMacro(vtkSQ3SubdivisionFilter, "$Revision: 1.20 $");
00030 vtkStandardNewMacro(vtkSQ3SubdivisionFilter);
00031 
00032 // Construct object with number of subdivisions set to 1.
00033 /*
00034 vtkSQ3SubdivisionFilter::vtkSQ3SubdivisionFilter()
00035 {
00036 }
00037 */
00038 
00039 void vtkSQ3SubdivisionFilter::Execute()
00040 {
00041   vtkIdType numPts, numCells;
00042   int level;
00043   vtkPoints *outputPts;
00044   vtkCellArray *outputPolys;
00045   vtkPolyData *input = this->GetInput();
00046   vtkPolyData *output = this->GetOutput();
00047   vtkPointData *outputPD;
00048   vtkCellData *outputCD;
00049 
00050   vtkDebugMacro(<< "Generating subdivision surface using interpolating scheme");
00051 
00052   
00053   if (input == NULL)
00054     {
00055     vtkErrorMacro(<<"Input is NULL");
00056     return;
00057     }
00058 
00059   numPts=input->GetNumberOfPoints();
00060   numCells=input->GetNumberOfCells();
00061 
00062   if (numPts < 1 || numCells < 1)
00063     {
00064     vtkDebugMacro(<<"No data to interpolate!");
00065     return;
00066     }
00067 
00068   //
00069   // Initialize and check input
00070   //
00071 
00072   vtkPolyData *inputDS = vtkPolyData::New();
00073   inputDS->CopyStructure (input);
00074   inputDS->GetPointData()->PassData(input->GetPointData());
00075   inputDS->GetCellData()->PassData(input->GetCellData());
00076 
00077   for (level = 0; level < this->NumberOfSubdivisions; level++)
00078     {
00079     // Generate topology  for the input dataset
00080     inputDS->BuildLinks();
00081     numCells = inputDS->GetNumberOfCells ();
00082 
00083     // Copy points from input. The new points will include the old points
00084     // and points calculated by the subdivision algorithm
00085     outputPts = vtkPoints::New();
00086     //outputPts->GetData()->DeepCopy(inputDS->GetPoints()->GetData());
00087                 outputPts->Allocate (numPts);
00088 
00089     // Copy pointdata structure from input
00090     outputPD = vtkPointData::New();
00091     outputPD->CopyAllocate(inputDS->GetPointData(),
00092                            2*inputDS->GetNumberOfPoints());
00093 
00094     // Copy celldata structure from input
00095     outputCD = vtkCellData::New();
00096     outputCD->CopyAllocate(inputDS->GetCellData(),3*numCells);
00097                 //outputCD->Print(*os);
00098 
00099     // Create triangles
00100     outputPolys = vtkCellArray::New();
00101     outputPolys->Allocate(outputPolys->EstimateSize(4*numCells,3));
00102 
00103                 // Do the split
00104     this->GenerateSubdivisionPoints (inputDS, outputPts, outputPD, outputPolys, outputCD);
00105 
00106     // start the next iteration with the input set to the output we just created
00107     inputDS->Delete();
00108     inputDS = vtkPolyData::New();
00109     inputDS->SetPoints(outputPts); outputPts->Delete();
00110     inputDS->SetPolys(outputPolys); outputPolys->Delete();
00111     inputDS->GetPointData()->PassData(outputPD); outputPD->Delete();
00112     inputDS->GetCellData()->PassData(outputCD); outputCD->Delete();
00113     inputDS->Squeeze();
00114     } // each level
00115 
00116   output->SetPoints(inputDS->GetPoints());
00117   output->SetPolys(inputDS->GetPolys());
00118   output->GetPointData()->PassData(inputDS->GetPointData());
00119   output->GetCellData()->PassData(inputDS->GetCellData());
00120   inputDS->Delete();
00121 }
00122 
00123 bool vtkSQ3SubdivisionFilter::generatePoint(vtkIdType cellId, vtkIdType &newId, vtkIdList *cellsVisited, 
00124                                                                                                                                                                                 vtkIdList *cellsPoints, vtkIdType *pts, vtkPolyData *inputDS, 
00125                                                                                                                                                                                 vtkPointData *outputPD, vtkPoints *outputPts)
00126 {
00127         //vtkIdType newId;
00128   vtkPoints *inputPts=inputDS->GetPoints();
00129   vtkIdList *pointIds = vtkIdList::New();
00130   vtkPointData *inputPD=inputDS->GetPointData();
00131         static float weights[3] = {.33333, .33333, .33334};
00132 
00133         pointIds->SetNumberOfIds(3);
00134         newId = cellsPoints->GetId(cellId);
00135 
00136         if ( cellsVisited->GetId(cellId) == 0 ) {
00137                 outputPD->CopyData(inputPD, pts[0], pts[0]);
00138                 outputPD->CopyData(inputPD, pts[1], pts[1]);
00139                 outputPD->CopyData(inputPD, pts[2], pts[2]);
00140 
00141                 pointIds->SetId(0,pts[0]);
00142                 pointIds->SetId(1,pts[1]);
00143                 pointIds->SetId(2,pts[2]);
00144                 newId = this->InterpolatePosition (inputPts, outputPts, pointIds, weights);
00145                 outputPD->InterpolatePoint (inputPD, newId, pointIds, weights);
00146 
00147                 cellsVisited->SetId(cellId, 1);
00148                 cellsPoints->SetId(cellId, newId);
00149 
00150                 return true;
00151 
00152                 /*
00153                 (*os) << "-------------------- Cell: " << cellId << "------------------\n";
00154                 //(*os) << ">>>> Edge Data:\n";
00155                 //edgeData->Print(*os);
00156                 (*os) << ">>>> Edge elements:\n";
00157                 for(int t=0; t<edgeData->GetNumberOfTuples(); ++t) {
00158                         for(int c=0; c<edgeData->GetNumberOfComponents(); ++c)
00159                                 (*os) << "Tuple[" << t << "]::Comp[" << c << "] = " << edgeData->GetComponent(t, c) << endl;
00160                 }
00161                 (*os) << ">>>> End of edge elements\n";
00162                 (*os) << ">>>> Points:\n";
00163                 outputPts->Print(*os);
00164                 for(int i=0; i<outputPts->GetNumberOfPoints(); ++i)
00165                         (*os) << "outputPts[" << i << "]= (" 
00166                                                 << outputPts->GetPoint(i)[0] << "," 
00167                                                 << outputPts->GetPoint(i)[1] << "," 
00168                                                 << outputPts->GetPoint(i)[2] << ")\n";
00169                 (*os) << ">>>> end of Points\n";
00170 
00171                 (*os) << endl << endl;
00172                 */
00173         } else {
00174                 return false;
00175         }
00176 }
00177 
00178 void vtkSQ3SubdivisionFilter::GenerateSubdivisionPoints (vtkPolyData *inputDS, 
00179                                                                                                                                                                                                                                  vtkPoints *outputPts, 
00180                                                                                                                                                                                                                                  vtkPointData *outputPD,
00181                                                                                                                                                                                                                                  vtkCellArray *outputPolys, 
00182                                                                                                                                                                                                                                  vtkCellData *outputCD)
00183 {
00184   vtkIdType                     *pts = 0;
00185   vtkIdType                     npts, cellId, newId;
00186   vtkCellArray  *inputPolys=inputDS->GetPolys();
00187   vtkIdList                     *cellIds = vtkIdList::New();
00188         vtkIdType                       cell = 0;
00189         vtkIdList                       *cellsVisited = vtkIdList::New();
00190         vtkIdList                       *cellsPoints = vtkIdList::New();
00191         vtkIdType                       newCellPts[3];
00192   vtkCellData           *inputCD = inputDS->GetCellData();
00193         vtkEdgeTable    *edgeTable = vtkEdgeTable::New();
00194         vtkIdType                       cellCount = inputDS->GetNumberOfCells();
00195         vtkIdList                       *stencil = vtkIdList::New();
00196         vtkPoints                       *inputPts=inputDS->GetPoints();
00197   vtkPointData  *inputPD=inputDS->GetPointData();
00198 
00199         cellsVisited->SetNumberOfIds(cellCount);
00200         cellsPoints->SetNumberOfIds(cellCount);
00201         for(int i=0; i<cellCount; ++i) {
00202                 cellsVisited->SetId(i, 0);
00203                 cellsPoints->SetId(i, -1);
00204         }
00205 
00206         vtkIdType numET = inputDS->GetNumberOfPoints() + inputDS->GetNumberOfCells();
00207         edgeTable->InitEdgeInsertion(numET);
00208 
00209         vtkIdType numPts = inputDS->GetNumberOfPoints();
00210         float *weights = new float[256];
00211         for (vtkIdType ptId=0; ptId < numPts; ptId++)
00212         {
00213                 this->GenerateEvenStencil (ptId, inputDS, stencil, weights);
00214     this->InterpolatePosition (inputPts, outputPts, stencil, weights);
00215                 outputPD->InterpolatePoint (inputPD, ptId, stencil, weights);
00216         }
00217 
00218   // Generate new points for subdivisions surface
00219   for (cellId=0, inputPolys->InitTraversal();
00220        inputPolys->GetNextCell(npts, pts); cellId++)
00221   {
00222     if ( inputDS->GetCellType(cellId) != VTK_TRIANGLE )
00223     {
00224       continue;
00225     }
00226 
00227                 bool actIsGen = generatePoint(cellId, newId, cellsVisited, cellsPoints, pts, inputDS, outputPD, outputPts);
00228 
00229                 //bool haveMoreCells = true;
00230                 vtkIdType p1 = pts[2], p2 = pts[0];
00231                 for (int edgeId=0; edgeId < 3; edgeId++) {
00232                         inputDS->GetCellEdgeNeighbors(cellId, p1, p2, cellIds);
00233                         if ( cellIds->GetNumberOfIds() != 0) {
00234                                 vtkIdType *neighpts = 0;
00235                                 vtkIdType neighnpts;
00236                                 vtkIdType newnId;
00237                                 inputDS->GetCellPoints(cellIds->GetId(cell), neighnpts, neighpts);
00238                                 bool neighIsGen = generatePoint(cellIds->GetId(cell), newnId, cellsVisited, cellsPoints, neighpts, inputDS, outputPD, outputPts);
00239 
00240                                 if ( edgeTable->IsEdge(newId, newnId) == -1 ) {         // Todo: Do a real edge existence check :)
00241                                         vtkIdType newCellPts[3];
00242                                         int id;
00243                                         vtkIdType myNewId;
00244 
00245                                         edgeTable->InsertEdge(newId, newnId);
00246 
00247                                         id = 0;
00248                                         newCellPts[id++] = p1;
00249                                         newCellPts[id++] = newnId;
00250                                         newCellPts[id++] = newId;
00251                                         myNewId = outputPolys->InsertNextCell (3, newCellPts);
00252                                         outputCD->CopyData(inputCD, cellId, newId);
00253                                         id = 0;
00254                                         newCellPts[id++] = newId;
00255                                         newCellPts[id++] = newnId;
00256                                         newCellPts[id++] = p2;
00257                                         myNewId = outputPolys->InsertNextCell (3, newCellPts);
00258                                         outputCD->CopyData(inputCD, cellId, newId);
00259                                 }
00260                         } else {
00261                                         vtkIdType newCellPts[3];
00262                                         int id;
00263                                         vtkIdType myNewId;
00264 
00265                                         id = 0;
00266                                         newCellPts[id++] = p1;
00267                                         newCellPts[id++] = p2;
00268                                         newCellPts[id++] = newId;
00269                                         myNewId = outputPolys->InsertNextCell (3, newCellPts);
00270                                         outputCD->CopyData(inputCD, cellId, myNewId);
00271                         }
00272 
00273                         p1 = p2;
00274                         if (edgeId < 2)
00275                         {
00276                                 p2 = pts[edgeId + 1];
00277                         }
00278                 }
00279 
00280   } // each cell
00281 
00282   cellIds->Delete();
00283         cellsVisited->Delete();
00284         cellsPoints->Delete();
00285         edgeTable->Delete();
00286 
00287 }
00288 
00289 void vtkSQ3SubdivisionFilter::GenerateSubdivisionCells (vtkPolyData *inputDS, 
00290                                                                                                                                                                                                                                 vtkIntArray *edgeData, 
00291                                                                                                                                                                                                                                 vtkCellArray *outputPolys, 
00292                                                                                                                                                                                                                                 vtkCellData *outputCD)
00293 {
00294   vtkIdType numCells = inputDS->GetNumberOfCells();
00295   vtkIdType cellId, newId;
00296   int id;
00297   vtkIdType npts;
00298   vtkIdType *pts;
00299   float edgePts[3];
00300   vtkIdType newCellPts[3];
00301   vtkCellData *inputCD = inputDS->GetCellData();
00302 
00303   // Now create new cells from existing points and generated edge points
00304   for (cellId=0; cellId < numCells; cellId++)
00305   {
00306     if ( inputDS->GetCellType(cellId) != VTK_TRIANGLE )
00307     {
00308       continue;
00309     }
00310     // get the original point ids and the ids stored as cell data
00311     inputDS->GetCellPoints(cellId, npts, pts);
00312     edgeData->GetTuple(cellId, edgePts);
00313     id = 0;
00314     newCellPts[id++] = pts[0];
00315     newCellPts[id++] = (int) edgePts[0];
00316     newCellPts[id++] = pts[2];
00317     newId = outputPolys->InsertNextCell (3, newCellPts);
00318     outputCD->CopyData (inputCD, cellId, newId);
00319 
00320     id = 0;
00321     newCellPts[id++] = pts[0];
00322     newCellPts[id++] = pts[1];
00323     newCellPts[id++] = (int) edgePts[0];
00324     newId = outputPolys->InsertNextCell (3, newCellPts);
00325     outputCD->CopyData (inputCD, cellId, newId);
00326 
00327     id = 0;
00328     newCellPts[id++] = pts[1];
00329     newCellPts[id++] = pts[2];
00330     newCellPts[id++] = (int) edgePts[0];
00331     newId = outputPolys->InsertNextCell (3, newCellPts);
00332     outputCD->CopyData (inputCD, cellId, newId);
00333         }
00334 }
00335 
00336 vtkIdType vtkSQ3SubdivisionFilter::InterpolatePosition (
00337         vtkPoints *inputPts, vtkPoints *outputPts,
00338         vtkIdList *stencil, float *weights)
00339 {
00340   float *xx, x[3];
00341   int i, j;
00342 
00343   for (j = 0; j < 3; j++)
00344     {
00345     x[j] = 0.0;
00346     }
00347 
00348   for (i = 0; i < stencil->GetNumberOfIds(); i++)
00349     {
00350     xx = inputPts->GetPoint(stencil->GetId(i));
00351     for (j = 0; j < 3; j++)
00352       {
00353       x[j] += xx[j] * weights[i];
00354       }
00355     }
00356   return outputPts->InsertNextPoint (x);
00357 
00358 }
00359 
00360 int vtkSQ3SubdivisionFilter::FindEdge (vtkPolyData *mesh,
00361                                                  vtkIdType cellId,
00362                                                  vtkIdType p1, vtkIdType p2,
00363                                                  vtkIntArray *edgeData,
00364                                                  vtkIdList *cellIds)
00365 {
00366   int edgeId = 0;
00367   int currentCellId = 0;
00368   int i;
00369   int numEdges;
00370   vtkIdType tp1, tp2;
00371   vtkCell *cell;
00372 
00373   // get all the cells that use the edge (except for cellId)
00374   mesh->GetCellEdgeNeighbors (cellId, p1, p2, cellIds);
00375 
00376   // find the edge that has the point we are looking for
00377   for ( i=0; i < cellIds->GetNumberOfIds(); i++)
00378     {
00379     currentCellId = cellIds->GetId(i);
00380     cell = mesh->GetCell(currentCellId);
00381     numEdges = cell->GetNumberOfEdges();
00382     tp1 = cell->GetPointId(2);
00383     tp2 = cell->GetPointId(0);
00384     for (edgeId=0; edgeId < numEdges; edgeId++)
00385       {
00386       if ( (tp1 == p1 && tp2 == p2) ||
00387            (tp2 == p1 && tp1 == p2))
00388         {
00389         break;
00390         }
00391       tp1 = tp2;
00392       tp2 = cell->GetPointId(edgeId + 1);
00393       }
00394     }
00395     // found the edge, return the stored value
00396     return (int) edgeData->GetComponent(currentCellId,edgeId);
00397 }
00398 
00399 void vtkSQ3SubdivisionFilter::GenerateEvenStencil (vtkIdType p1,
00400                                                     vtkPolyData *polys,
00401                                                     vtkIdList *stencilIds,
00402                                                     float *weights)
00403 {
00404   vtkIdList *cellIds = vtkIdList::New();
00405   vtkIdList *ptIds = vtkIdList::New();
00406   vtkCell *cell;
00407 
00408   int i, j;
00409   int numCellsInLoop;
00410   int startCell, nextCell;
00411   vtkIdType p, p2;
00412   vtkIdType bp1, bp2;
00413   int K;
00414   float beta, cosSQ;
00415 
00416   // Get the cells that use this point
00417   polys->GetPointCells (p1, cellIds);
00418   numCellsInLoop = cellIds->GetNumberOfIds();
00419   if (numCellsInLoop < 1)
00420       {
00421       vtkWarningMacro("numCellsInLoop < 1: " << numCellsInLoop);
00422       stencilIds->Reset();
00423       return;
00424       }
00425   // Find an edge to start with that contains p1
00426   polys->GetCellPoints (cellIds->GetId(0), ptIds);
00427   p2 = ptIds->GetId(0);
00428   i = 1;
00429   while (p1 == p2)
00430     {
00431     p2 = ptIds->GetId(i++);
00432     }
00433   polys->GetCellEdgeNeighbors (-1, p1, p2, cellIds);
00434 
00435   nextCell = cellIds->GetId(0);
00436   bp2 = -1;
00437   bp1 = p2;
00438   if (cellIds->GetNumberOfIds() == 1)
00439     {
00440     startCell = -1;
00441     }
00442   else
00443     {
00444     startCell = cellIds->GetId(1);
00445     }
00446 
00447   stencilIds->Reset();
00448   stencilIds->InsertNextId(p2);
00449 
00450   // walk around the loop counter-clockwise and get cells
00451   for (j = 0; j < numCellsInLoop; j++)
00452     {
00453     cell = polys->GetCell(nextCell);
00454     p = -1;
00455     for (i = 0; i < 3; i++)
00456       {
00457       if ((p = cell->GetPointId(i)) != p1 && cell->GetPointId(i) != p2)
00458         {
00459         break;
00460         }
00461       }
00462     p2 = p;
00463     stencilIds->InsertNextId (p2);
00464     polys->GetCellEdgeNeighbors (nextCell, p1, p2, cellIds);
00465     if (cellIds->GetNumberOfIds() != 1)
00466       {
00467       bp2 = p2;
00468       j++;
00469       break;
00470       }
00471     nextCell = cellIds->GetId(0);
00472     }
00473 
00474   // now walk around the other way. this will only happen if there
00475   // is a boundary cell left that we have not visited
00476   nextCell = startCell;
00477   p2 = bp1;
00478   for (; j < numCellsInLoop && startCell != -1; j++)
00479     {
00480     cell = polys->GetCell(nextCell);
00481     p = -1;
00482     for (i = 0; i < 3; i++)
00483       {
00484       if ((p = cell->GetPointId(i)) != p1 && cell->GetPointId(i) != p2)
00485         {
00486         break;
00487         }
00488       }
00489     p2 = p;
00490     stencilIds->InsertNextId (p2);
00491     polys->GetCellEdgeNeighbors (nextCell, p1, p2, cellIds);
00492     if (cellIds->GetNumberOfIds() != 1)
00493       {
00494       bp1 = p2;
00495       break;
00496       }
00497     nextCell = cellIds->GetId(0);
00498     }
00499 
00500   if (bp2 != -1) // boundary edge
00501     {
00502     stencilIds->SetNumberOfIds(3);
00503     stencilIds->SetId(0,bp2);
00504     stencilIds->SetId(1,bp1);
00505     stencilIds->SetId(2,p1);
00506     weights[0] = .125;
00507     weights[1] = .125;
00508     weights[2] = .75;
00509     }
00510   else
00511     {
00512     K = stencilIds->GetNumberOfIds();
00513    // Remove last id. It's a duplicate of the first
00514     K--;
00515     if (K > 3)
00516       {
00517       // Generate weights
00518 #define VTK_PI 3.14159265358979
00519       cosSQ = .375 + .25 * cos (2.0 * VTK_PI / (float) K);
00520       cosSQ = cosSQ * cosSQ;
00521       beta = (.625 -  cosSQ) / (float) K;
00522       }
00523     else
00524       {
00525       beta = 3.0 / 16.0;
00526       }
00527     for (j = 0; j < K; j++)
00528       {
00529       weights[j] = beta;
00530       }
00531     weights[K] = 1.0 - K * beta;
00532     stencilIds->SetId (K,p1);
00533     }
00534   cellIds->Delete();
00535   ptIds->Delete();
00536 }
00537 
00538 void vtkSQ3SubdivisionFilter::PrintSelf(ostream& os, vtkIndent indent)
00539 {
00540   this->Superclass::PrintSelf(os,indent);
00541 
00542   os << indent << "Number of subdivisions: " << this->NumberOfSubdivisions << endl;
00543 }
00544 
00545