/*========================================================================= Program: Visualization Toolkit Module: vtkQuadricClustering.cxx Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen All rights reserved. See Copyright.txt or http://www.kitware.com/Copyright.htm for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notice for more information. =========================================================================*/ #include "vtkQuadricClustering.h" #include "vtkCellArray.h" #include "vtkCellData.h" #include "vtkExecutive.h" #include "vtkFeatureEdges.h" #include "vtkInformation.h" #include "vtkInformationVector.h" #include "vtkMath.h" #include "vtkObjectFactory.h" #include "vtkPointData.h" #include "vtkPolyData.h" #include "vtkTimerLog.h" #include "vtkTriangle.h" #include // keep track of inserted triangles vtkStandardNewMacro(vtkQuadricClustering); //---------------------------------------------------------------------------- // PIMPLd STL set for keeping track of inserted cells struct vtkQuadricClusteringIdTypeHash { size_t operator()(vtkIdType val) const { return static_cast(val); } }; class vtkQuadricClusteringCellSet : public vtksys::hash_set {}; typedef vtkQuadricClusteringCellSet::iterator vtkQuadricClusteringCellSetIterator; //---------------------------------------------------------------------------- // Construct with default NumberOfDivisions to 50, DivisionSpacing to 1 // in all (x,y,z) directions. AutoAdjustNumberOfDivisions is set to ON. // ComputeNumberOfDivisions to OFF. UseFeatureEdges and UseFeaturePoints // are set to OFF by default // The default behavior is also to compute an optimal position in each // bin to produce the output triangles (this is also recommended) vtkQuadricClustering::vtkQuadricClustering() { this->Bounds[0] = this->Bounds[1] = this->Bounds[2] = 0.0; this->Bounds[3] = this->Bounds[4] = this->Bounds[5] = 0.0; this->NumberOfXDivisions = 50; this->NumberOfYDivisions = 50; this->NumberOfZDivisions = 50; this->QuadricArray = NULL; this->NumberOfBinsUsed = 0; this->AbortExecute = 0; this->AutoAdjustNumberOfDivisions = 1; this->ComputeNumberOfDivisions = 0; this->DivisionOrigin[0] = 0.0; this->DivisionOrigin[1] = 0.0; this->DivisionOrigin[2] = 0.0; this->DivisionSpacing[0] = 1.0; this->DivisionSpacing[1] = 1.0; this->DivisionSpacing[2] = 1.0; this->UseFeatureEdges = 0; this->UseFeaturePoints = 0; this->FeaturePointsAngle = 30.0; this->UseInternalTriangles = 1; this->UseInputPoints = 0; this->PreventDuplicateCells = 1; this->CellSet = NULL; this->NumberOfBins = 0; this->OutputTriangleArray = NULL; this->OutputLines = NULL; // Used for matching boundaries. this->FeatureEdges = vtkFeatureEdges::New(); this->FeatureEdges->FeatureEdgesOff(); this->FeatureEdges->BoundaryEdgesOn(); this->FeaturePoints = vtkPoints::New(); this->InCellCount = this->OutCellCount = 0; this->CopyCellData = 0; this->GetInformation()->Set(vtkAlgorithm::PRESERVES_RANGES(), 1); this->GetInformation()->Set(vtkAlgorithm::PRESERVES_BOUNDS(), 1); } //---------------------------------------------------------------------------- vtkQuadricClustering::~vtkQuadricClustering() { this->FeatureEdges->Delete(); this->FeatureEdges = NULL; this->FeaturePoints->Delete(); this->FeaturePoints = NULL; if (this->CellSet) { delete this->CellSet; this->CellSet = NULL; } if (this->QuadricArray) { delete [] this->QuadricArray; this->QuadricArray = NULL; } if (this->OutputTriangleArray) { this->OutputTriangleArray->Delete(); this->OutputTriangleArray = NULL; } if (this->OutputLines) { this->OutputLines->Delete(); this->OutputLines = NULL; } } //---------------------------------------------------------------------------- int vtkQuadricClustering::RequestData( vtkInformation *vtkNotUsed(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 vtkPolyData *input = 0; if (inInfo) { input = vtkPolyData::SafeDownCast( inInfo->Get(vtkDataObject::DATA_OBJECT())); } vtkPolyData *output = vtkPolyData::SafeDownCast( outInfo->Get(vtkDataObject::DATA_OBJECT())); vtkTimerLog *tlog=NULL; if (!input || (input->GetNumberOfPoints() == 0)) { // The user may be calling StartAppend, Append, and EndAppend explicitly. return 1; } if (input->CheckAttributes()) { // avoid crashing if input is not all we expect (is not consistent). return 1; } if (this->Debug) { tlog = vtkTimerLog::New(); tlog->StartTimer(); } // Lets limit the number of divisions based on // the number of points in the input. int target = input->GetNumberOfPoints(); int numDiv = (this->NumberOfXDivisions * this->NumberOfYDivisions * this->NumberOfZDivisions) / 2; if (this->AutoAdjustNumberOfDivisions && numDiv > target) { double factor = pow(((double)numDiv/(double)target),0.33333); this->NumberOfDivisions[0] = (int)(0.5+(double)(this->NumberOfXDivisions)/factor); this->NumberOfDivisions[0] = (this->NumberOfDivisions[0] > 0 ? this->NumberOfDivisions[0] : 1); this->NumberOfDivisions[1] = (int)(0.5+(double)(this->NumberOfYDivisions)/factor); this->NumberOfDivisions[1] = (this->NumberOfDivisions[1] > 0 ? this->NumberOfDivisions[1] : 1); this->NumberOfDivisions[2] = (int)(0.5+(double)(this->NumberOfZDivisions)/factor); this->NumberOfDivisions[2] = (this->NumberOfDivisions[2] > 0 ? this->NumberOfDivisions[2] : 1); } else { this->NumberOfDivisions[0] = this->NumberOfXDivisions; this->NumberOfDivisions[1] = this->NumberOfYDivisions; this->NumberOfDivisions[2] = this->NumberOfZDivisions; } this->UpdateProgress(.01); this->StartAppend(input->GetBounds()); this->UpdateProgress(.2); this->SliceSize = this->NumberOfDivisions[0]*this->NumberOfDivisions[1]; this->Append(input); if (this->UseFeatureEdges) { // Adjust bin points that contain boundary edges. this->AppendFeatureQuadrics(input, output); } if (this->UseInputPoints) { this->EndAppendUsingPoints(input, output); } else { this->EndAppend(); } // Free up some memory. if (this->QuadricArray) { delete [] this->QuadricArray; this->QuadricArray = NULL; } if ( this->Debug ) { tlog->StopTimer(); vtkDebugMacro(<<"Execution took: "<GetElapsedTime()<<" seconds."); tlog->Delete(); } return 1; } //---------------------------------------------------------------------------- void vtkQuadricClustering::StartAppend(double *bounds) { vtkIdType i; // If there are duplicate triangles. remove them if ( this->PreventDuplicateCells ) { this->CellSet = new vtkQuadricClusteringCellSet; this->NumberOfBins = this->NumberOfDivisions[0]*this->NumberOfDivisions[1]*this->NumberOfDivisions[2]; } // Copy over the bounds. for (i = 0; i < 6; ++i) { this->Bounds[i]= bounds[i]; } if (this->ComputeNumberOfDivisions) { // extend the bounds so that it will not produce fractions of bins. double x, y, z; x = floor((bounds[0]-this->DivisionOrigin[0])/this->DivisionSpacing[0]); y = floor((bounds[2]-this->DivisionOrigin[1])/this->DivisionSpacing[1]); z = floor((bounds[4]-this->DivisionOrigin[2])/this->DivisionSpacing[2]); this->Bounds[0] = this->DivisionOrigin[0]+(x * this->DivisionSpacing[0]); this->Bounds[2] = this->DivisionOrigin[1]+(y * this->DivisionSpacing[1]); this->Bounds[4] = this->DivisionOrigin[2]+(z * this->DivisionSpacing[2]); x = ceil((bounds[1]-this->Bounds[0])/this->DivisionSpacing[0]); y = ceil((bounds[3]-this->Bounds[2])/this->DivisionSpacing[1]); z = ceil((bounds[5]-this->Bounds[4])/this->DivisionSpacing[2]); this->Bounds[1] = this->Bounds[0] + (x * this->DivisionSpacing[0]); this->Bounds[3] = this->Bounds[2] + (y * this->DivisionSpacing[1]); this->Bounds[5] = this->Bounds[4] + (z * this->DivisionSpacing[2]); this->NumberOfDivisions[0] = (int)x > 0 ? (int)x : 1; this->NumberOfDivisions[1] = (int)y > 0 ? (int)y : 1; this->NumberOfDivisions[2] = (int)z > 0 ? (int)z : 1; } else { this->DivisionOrigin[0] = bounds[0]; this->DivisionOrigin[1] = bounds[2]; this->DivisionOrigin[2] = bounds[4]; this->DivisionSpacing[0] = (bounds[1]-bounds[0])/this->NumberOfDivisions[0]; this->DivisionSpacing[1] = (bounds[3]-bounds[2])/this->NumberOfDivisions[1]; this->DivisionSpacing[2] = (bounds[5]-bounds[4])/this->NumberOfDivisions[2]; } // Check for conditions that can occur if the Append methods // are not called in the correct order. if (this->OutputTriangleArray) { this->OutputTriangleArray->Delete(); this->OutputTriangleArray = NULL; //vtkWarningMacro("Array already created. Did you call EndAppend?"); } if (this->OutputLines) { this->OutputLines->Delete(); this->OutputLines = NULL; //vtkWarningMacro("Array already created. Did you call EndAppend?"); } this->OutputTriangleArray = vtkCellArray::New(); this->OutputLines = vtkCellArray::New(); this->XBinSize = (this->Bounds[1]-this->Bounds[0])/this->NumberOfDivisions[0]; this->YBinSize = (this->Bounds[3]-this->Bounds[2])/this->NumberOfDivisions[1]; this->ZBinSize = (this->Bounds[5]-this->Bounds[4])/this->NumberOfDivisions[2]; this->XBinStep = (this->XBinSize > 0.0) ? (1.0/this->XBinSize) : 0.0; this->YBinStep = (this->YBinSize > 0.0) ? (1.0/this->YBinSize) : 0.0; this->ZBinStep = (this->ZBinSize > 0.0) ? (1.0/this->ZBinSize) : 0.0; this->NumberOfBinsUsed = 0; if (this->QuadricArray) { delete [] this->QuadricArray; this->QuadricArray = NULL; } this->QuadricArray = new vtkQuadricClustering::PointQuadric[this->NumberOfDivisions[0] * this->NumberOfDivisions[1] * this->NumberOfDivisions[2]]; if (this->QuadricArray == NULL) { vtkErrorMacro("Could not allocate quadric grid."); return; } vtkInformation *inInfo = this->GetExecutive()->GetInputInformation(0, 0); vtkInformation *outInfo = this->GetExecutive()->GetOutputInformation(0); vtkPolyData *input = 0; if (inInfo) { input = vtkPolyData::SafeDownCast( inInfo->Get(vtkDataObject::DATA_OBJECT())); } vtkPolyData *output = vtkPolyData::SafeDownCast( outInfo->Get(vtkDataObject::DATA_OBJECT())); // Allocate CellData here. if (this->CopyCellData && input) { output->GetCellData()->CopyAllocate( input->GetCellData(), this->NumberOfBinsUsed); this->InCellCount = this->OutCellCount = 0; } } //---------------------------------------------------------------------------- void vtkQuadricClustering::Append(vtkPolyData *pd) { vtkCellArray *inputPolys, *inputStrips, *inputLines, *inputVerts; vtkPoints *inputPoints = pd->GetPoints(); // Check for mis-use of the Append methods. if (this->OutputTriangleArray == NULL || this->OutputLines == NULL) { vtkDebugMacro("Missing Array: Did you call StartAppend?"); return; } vtkInformation *outInfo = this->GetExecutive()->GetOutputInformation(0); vtkPolyData *output = vtkPolyData::SafeDownCast( outInfo->Get(vtkDataObject::DATA_OBJECT())); inputVerts = pd->GetVerts(); if (inputVerts) { this->AddVertices(inputVerts, inputPoints, 1, pd, output); } this->UpdateProgress(.40); inputLines = pd->GetLines(); if (inputLines) { this->AddEdges(inputLines, inputPoints, 1, pd, output); } this->UpdateProgress(.60); inputPolys = pd->GetPolys(); if (inputPolys) { this->AddPolygons(inputPolys, inputPoints, 1, pd, output); } this->UpdateProgress(.80); inputStrips = pd->GetStrips(); if (inputStrips) { this->AddStrips(inputStrips, inputPoints, 1, pd, output); } } //---------------------------------------------------------------------------- void vtkQuadricClustering::AddPolygons(vtkCellArray *polys, vtkPoints *points, int geometryFlag, vtkPolyData *input, vtkPolyData *output) { int j; vtkIdType *ptIds = 0; vtkIdType numPts = 0; double pts0[3], pts1[3], pts2[3]; vtkIdType binIds[3]; double total = polys->GetNumberOfCells(); double curr = 0; double step = total / 10; if (step < 1000.0) { step = 1000.0; } double cstep = step; for ( polys->InitTraversal(); polys->GetNextCell(numPts, ptIds); ) { points->GetPoint(ptIds[0],pts0); binIds[0] = this->HashPoint(pts0); for (j=0; j < numPts-2; j++)//creates triangles; assumes poly is convex { points->GetPoint(ptIds[j+1],pts1); binIds[1] = this->HashPoint(pts1); points->GetPoint(ptIds[j+2],pts2); binIds[2] = this->HashPoint(pts2); this->AddTriangle(binIds, pts0, pts1, pts2, geometryFlag, input, output); } ++this->InCellCount; if ( curr > cstep ) { this->UpdateProgress(.6 + .2 * curr / total); cstep += step; } curr += 1; }//for all polygons } //---------------------------------------------------------------------------- void vtkQuadricClustering::AddStrips(vtkCellArray *strips, vtkPoints *points, int geometryFlag, vtkPolyData *input, vtkPolyData *output) { int j; vtkIdType *ptIds = 0; vtkIdType numPts = 0; double pts[3][3]; vtkIdType binIds[3]; int odd; // Used to flip order of every other triangle in a strip. for ( strips->InitTraversal(); strips->GetNextCell(numPts, ptIds); ) { points->GetPoint(ptIds[0], pts[0]); binIds[0] = this->HashPoint(pts[0]); points->GetPoint(ptIds[1], pts[1]); binIds[1] = this->HashPoint(pts[1]); // This internal loop handles triangle strips. odd = 0; for (j = 2; j < numPts; ++j) { points->GetPoint(ptIds[j], pts[2]); binIds[2] = this->HashPoint(pts[2]); this->AddTriangle(binIds, pts[0], pts[1], pts[2], geometryFlag, input, output); pts[odd][0] = pts[2][0]; pts[odd][1] = pts[2][1]; pts[odd][2] = pts[2][2]; binIds[odd] = binIds[2]; // Toggle odd. odd = odd ? 0 : 1; } ++this->InCellCount; } } //---------------------------------------------------------------------------- inline void vtkQuadricClustering::InitializeQuadric(double quadric[9]) { quadric[0] = 0.0; quadric[1] = 0.0; quadric[2] = 0.0; quadric[3] = 0.0; quadric[4] = 0.0; quadric[5] = 0.0; quadric[6] = 0.0; quadric[7] = 0.0; quadric[8] = 0.0; } //---------------------------------------------------------------------------- // The error function is the volume (squared) of the tetrahedron formed by the // triangle and the point. We ignore constant factors across all coefficents, // and the constant coefficient. // If geomertyFlag is 1 then the triangle is added to the output. Otherwise, // only the quadric is affected. void vtkQuadricClustering::AddTriangle(vtkIdType *binIds, double *pt0, double *pt1, double *pt2, int geometryFlag, vtkPolyData *input, vtkPolyData *output) { int i; vtkIdType triPtIds[3]; double quadric[9], quadric4x4[4][4]; vtkIdType minIdx, midIdx, maxIdx, idx; // Special condition for fast execution. // Only add triangles that traverse three bins to quadrics. if (this->UseInternalTriangles == 0) { if (binIds[0] == binIds[1] || binIds[0] == binIds[2] || binIds[1] == binIds[2]) { return; } } // Compute the quadric. vtkTriangle::ComputeQuadric(pt0, pt1, pt2, quadric4x4); quadric[0] = quadric4x4[0][0]; quadric[1] = quadric4x4[0][1]; quadric[2] = quadric4x4[0][2]; quadric[3] = quadric4x4[0][3]; quadric[4] = quadric4x4[1][1]; quadric[5] = quadric4x4[1][2]; quadric[6] = quadric4x4[1][3]; quadric[7] = quadric4x4[2][2]; quadric[8] = quadric4x4[2][3]; // Add the quadric to each of the three corner bins. for (i = 0; i < 3; ++i) { // If the current quadric is not initialized, then clear it out. if (this->QuadricArray[binIds[i]].Dimension > 2) { this->QuadricArray[binIds[i]].Dimension = 2; // Initialize the coeff this->InitializeQuadric(this->QuadricArray[binIds[i]].Quadric); } if (this->QuadricArray[binIds[i]].Dimension == 2) { // Points and segments supercede triangles. this->AddQuadric(binIds[i], quadric); } } if (geometryFlag) { // Now add the triangle to the geometry. for (i = 0; i < 3; i++) { // Get the vertex from each bin. if (this->QuadricArray[binIds[i]].VertexId == -1) { this->QuadricArray[binIds[i]].VertexId = this->NumberOfBinsUsed; this->NumberOfBinsUsed++; } triPtIds[i] = this->QuadricArray[binIds[i]].VertexId; } // This comparison could just as well be on triPtIds. if (binIds[0] != binIds[1] && binIds[0] != binIds[2] && binIds[1] != binIds[2]) { if ( this->PreventDuplicateCells ) { minIdx = ( binIds[0] binIds[2] ) { maxIdx = 1; midIdx = 2; } else { maxIdx = 2; midIdx = 1; } break; case 1: if ( binIds[0] > binIds[2] ) { maxIdx = 0; midIdx = 2; } else { maxIdx = 2; midIdx = 0; } break; case 2: if ( binIds[0] > binIds[1] ) { maxIdx = 0; midIdx = 1; } else { maxIdx = 1; midIdx = 0; } break; } idx = binIds[minIdx] + this->NumberOfBins*binIds[midIdx] + this->NumberOfBins*this->NumberOfBins*binIds[maxIdx]; if ( this->CellSet->find(idx) == this->CellSet->end() ) { this->CellSet->insert(idx); this->OutputTriangleArray->InsertNextCell(3, triPtIds); if (this->CopyCellData && input) { output->GetCellData()-> CopyData(input->GetCellData(), this->InCellCount,this->OutCellCount++); }//if cell data }//if not a duplicate } else //don't check for duplicates { this->OutputTriangleArray->InsertNextCell(3, triPtIds); if (this->CopyCellData && input) { output->GetCellData()-> CopyData(input->GetCellData(), this->InCellCount,this->OutCellCount++); }//if cell data }//don't check for duplicates }//if not duplicate vertices }//if this should be inserted } //---------------------------------------------------------------------------- void vtkQuadricClustering::AddEdges(vtkCellArray *edges, vtkPoints *points, int geometryFlag, vtkPolyData *input, vtkPolyData *output) { int j; vtkIdType numCells, i; vtkIdType *ptIds = 0; vtkIdType numPts = 0; double pt0[3], pt1[3]; vtkIdType binIds[2]; // Add the edges to the error fuction. numCells = edges->GetNumberOfCells(); edges->InitTraversal(); for (i = 0; i < numCells; ++i) { edges->GetNextCell(numPts, ptIds); if(numPts != 0) { points->GetPoint(ptIds[0], pt0); binIds[0] = this->HashPoint(pt0); // This internal loop handles line strips. for (j = 1; j < numPts; ++j) { points->GetPoint(ptIds[j], pt1); binIds[1] = this->HashPoint(pt1); this->AddEdge(binIds, pt0, pt1, geometryFlag, input, output); pt0[0] = pt1[0]; pt0[1] = pt1[1]; pt0[2] = pt1[2]; binIds[0] = binIds[1]; } } ++this->InCellCount; } } //---------------------------------------------------------------------------- // The error function is the square of the area of the triangle formed by the // edge and the point. We ignore constants across all terms. // If geometryFlag is 1 then the edge is added to the output. Otherwise, // only the quadric is affected. void vtkQuadricClustering::AddEdge(vtkIdType *binIds, double *pt0, double *pt1, int geometryFlag, vtkPolyData *input, vtkPolyData *output) { int i; vtkIdType edgePtIds[2]; double length2, tmp; double d[3]; double m[3]; // The mid point of the segement.(p1 or p2 could be used also). double md; // The dot product of m and d. double q[9]; // Compute quadric for line segment. // Line segment quadric is the area (squared) of the triangle (seg,pt) // Compute the direction vector of the segment. d[0] = pt1[0] - pt0[0]; d[1] = pt1[1] - pt0[1]; d[2] = pt1[2] - pt0[2]; // Compute the length^2 of the line segement. length2 = d[0]*d[0] + d[1]*d[1] + d[2]*d[2]; if (length2 == 0.0) { // Coincident points. Avoid divide by zero. return; } // Normalize the direction vector. tmp = 1.0 / sqrt(length2); d[0] = d[0] * tmp; d[1] = d[1] * tmp; d[2] = d[2] * tmp; // Compute the mid point of the segment. m[0] = 0.5 * (pt1[0] + pt0[0]); m[1] = 0.5 * (pt1[1] + pt0[1]); m[2] = 0.5 * (pt1[2] + pt0[2]); // Compute dot(m, d); md = m[0]*d[0] + m[1]*d[1] + m[2]*d[2]; // We save nine coefficients of the error function cooresponding to: // 0: Px^2 // 1: PxPy // 2: PxPz // 3: Px // 4: Py^2 // 5: PyPz // 6: Py // 7: Pz^2 // 8: Pz // We ignore the constant because it disappears with the derivative. q[0] = length2*(1.0 - d[0]*d[0]); q[1] = -length2*(d[0]*d[1]); q[2] = -length2*(d[0]*d[2]); q[3] = length2*(d[0]*md - m[0]); q[4] = length2*(1.0 - d[1]*d[1]); q[5] = -length2*(d[1]*d[2]); q[6] = length2*(d[1]*md - m[1]); q[7] = length2*(1.0 - d[2]*d[2]); q[8] = length2*(d[2]*md - m[2]); for (i = 0; i < 2; ++i) { // If the current quadric is from triangles (or not initialized), then clear it out. if (this->QuadricArray[binIds[i]].Dimension > 1) { this->QuadricArray[binIds[i]].Dimension = 1; // Initialize the coeff this->InitializeQuadric(this->QuadricArray[binIds[i]].Quadric); } if (this->QuadricArray[binIds[i]].Dimension == 1) { // Points supercede segements. this->AddQuadric(binIds[i], q); } } if (geometryFlag) { // Now add the edge to the geometry. for (i = 0; i < 2; i++) { // Get the vertex from each bin. if (this->QuadricArray[binIds[i]].VertexId == -1) { this->QuadricArray[binIds[i]].VertexId = this->NumberOfBinsUsed; this->NumberOfBinsUsed++; } edgePtIds[i] = this->QuadricArray[binIds[i]].VertexId; } // This comparison could just as well be on edgePtIds. if (binIds[0] != binIds[1]) { this->OutputLines->InsertNextCell(2, edgePtIds); if (this->CopyCellData && input) { output->GetCellData()-> CopyData(input->GetCellData(),this->InCellCount, this->OutCellCount++); } } } } //---------------------------------------------------------------------------- void vtkQuadricClustering::AddVertices(vtkCellArray *verts, vtkPoints *points, int geometryFlag, vtkPolyData *input, vtkPolyData *output) { int j; vtkIdType numCells, i; vtkIdType *ptIds = 0; vtkIdType numPts = 0; double pt[3]; vtkIdType binId; numCells = verts->GetNumberOfCells(); double cstep = (double)numCells / 10.0; if (cstep < 1000.0) { cstep = 1000.0; } double next = cstep; double curr = 0; verts->InitTraversal(); for (i = 0; i < numCells; ++i) { verts->GetNextCell(numPts, ptIds); // Can there be poly vertices? for (j = 0; j < numPts; ++j) { points->GetPoint(ptIds[j], pt); binId = this->HashPoint(pt); this->AddVertex(binId, pt, geometryFlag, input, output); } ++this->InCellCount; if ( curr > next ) { this->UpdateProgress(.2 + .2 * curr / (double)numCells); next += cstep; } curr += 1; } } //---------------------------------------------------------------------------- // The error function is the length (point to vert) squared. // We ignore constants across all terms. // If geomertyFlag is 1 then the vert is added to the output. Otherwise, // only the quadric is affected. void vtkQuadricClustering::AddVertex(vtkIdType binId, double *pt, int geometryFlag, vtkPolyData *input, vtkPolyData *output) { double q[9]; // Compute quadric for the vertex. // We save nine coefficients of the error function cooresponding to: // 0: Px^2 // 1: PxPy // 2: PxPz // 3: Px // 4: Py^2 // 5: PyPz // 6: Py // 7: Pz^2 // 8: Pz // We ignore the constant because it disappears with the derivative. q[0] = 1.0; q[1] = 0.0; q[2] = 0.0; q[3] = -pt[0]; q[4] = 1.0; q[5] = 0.0; q[6] = -pt[1]; q[7] = 1.0; q[8] = -pt[2]; // If the current quadric is from triangles, edges (or not initialized), // then clear it out. if (this->QuadricArray[binId].Dimension > 0) { this->QuadricArray[binId].Dimension = 0; // Initialize the coeff this->InitializeQuadric(this->QuadricArray[binId].Quadric); } if (this->QuadricArray[binId].Dimension == 0) { // Points supercede all other types of quadrics. this->AddQuadric(binId, q); } if (geometryFlag) { // Now add the vert to the geometry. // Get the vertex from the bin. if (this->QuadricArray[binId].VertexId == -1) { this->QuadricArray[binId].VertexId = this->NumberOfBinsUsed; this->NumberOfBinsUsed++; if (this->CopyCellData && input) { output->GetCellData()-> CopyData(input->GetCellData(), this->InCellCount, this->OutCellCount++); } } } } //---------------------------------------------------------------------------- void vtkQuadricClustering::AddQuadric(vtkIdType binId, double quadric[9]) { double *q = this->QuadricArray[binId].Quadric; for (int i=0; i<9; i++) { q[i] += (quadric[i] * 100000000.0); } } //---------------------------------------------------------------------------- vtkIdType vtkQuadricClustering::HashPoint(double point[3]) { vtkIdType binId; int xBinCoord; int yBinCoord; int zBinCoord; xBinCoord = static_cast((point[0] - this->Bounds[0])*this->XBinStep); if (xBinCoord < 0) { xBinCoord = 0; } else if (xBinCoord >= this->NumberOfDivisions[0]) { xBinCoord = this->NumberOfDivisions[0] - 1; } yBinCoord = static_cast((point[1] - this->Bounds[2])*this->YBinStep); if (yBinCoord < 0) { yBinCoord = 0; } else if (yBinCoord >= this->NumberOfDivisions[1]) { yBinCoord = this->NumberOfDivisions[1] - 1; } zBinCoord = static_cast((point[2] - this->Bounds[4])*this->ZBinStep); if (zBinCoord < 0) { zBinCoord = 0; } else if (zBinCoord >= this->NumberOfDivisions[2]) { zBinCoord = this->NumberOfDivisions[2] - 1; } // vary x fastest, then y, then z binId = xBinCoord + yBinCoord*this->NumberOfDivisions[0] + zBinCoord*this->SliceSize; return binId; } //---------------------------------------------------------------------------- void vtkQuadricClustering::EndAppend() { vtkInformation *inInfo = this->GetExecutive()->GetInputInformation(0, 0); vtkInformation *outInfo = this->GetExecutive()->GetOutputInformation(0); vtkPolyData *input = 0; if (inInfo) { input = vtkPolyData::SafeDownCast( inInfo->Get(vtkDataObject::DATA_OBJECT())); } vtkPolyData *output = vtkPolyData::SafeDownCast( outInfo->Get(vtkDataObject::DATA_OBJECT())); vtkIdType i, numBuckets; int abortExecute=0; vtkPoints *outputPoints; double newPt[3]; numBuckets = this->NumberOfDivisions[0] * this->NumberOfDivisions[1] * this->NumberOfDivisions[2]; double step = (double)numBuckets / 10.0; if (step < 1000.0) { step = 1000.0; } double cstep = 0; // Check for mis use of the Append methods. if (this->OutputTriangleArray == NULL || this->OutputLines == NULL) { vtkDebugMacro("Missing Array: Did you call StartAppend?"); return; } // Clean up if ( this->PreventDuplicateCells ) { delete this->CellSet; this->CellSet = NULL; } // Compute the representative points for each bin outputPoints = vtkPoints::New(); for (i = 0; !abortExecute && i < numBuckets; i++ ) { if (cstep > step) { cstep = 0; vtkDebugMacro(<<"Finding point in bin #" << i); this->UpdateProgress (0.8+0.2*i/numBuckets); abortExecute = this->GetAbortExecute(); } ++cstep; if (this->QuadricArray[i].VertexId != -1) { this->ComputeRepresentativePoint(this->QuadricArray[i].Quadric, i, newPt); outputPoints->InsertPoint(this->QuadricArray[i].VertexId, newPt); } } // Set up the output data object. output->SetPoints(outputPoints); outputPoints->Delete(); if (this->OutputTriangleArray->GetNumberOfCells() > 0) { output->SetPolys(this->OutputTriangleArray); } this->OutputTriangleArray->Delete(); this->OutputTriangleArray = NULL; if (this->OutputLines->GetNumberOfCells() > 0) { output->SetLines(this->OutputLines); } this->OutputLines->Delete(); this->OutputLines = NULL; this->EndAppendVertexGeometry(input, output); // Tell the data is is up to date // (in case the user calls this method directly). output->DataHasBeenGenerated(); // Free the quadric array. if (this->QuadricArray) { delete [] this->QuadricArray; this->QuadricArray = NULL; } } //---------------------------------------------------------------------------- void vtkQuadricClustering::ComputeRepresentativePoint(double quadric[9], vtkIdType binId, double point[3]) { int i, j; double A[3][3], U[3][3], UT[3][3], VT[3][3], V[3][3]; double b[3], w[3]; double W[3][3], tempMatrix[3][3]; double cellCenter[3], tempVector[3]; double cellBounds[6]; int x, y, z; double quadric4x4[4][4]; quadric4x4[0][0] = quadric[0]; quadric4x4[0][1] = quadric4x4[1][0] = quadric[1]; quadric4x4[0][2] = quadric4x4[2][0] = quadric[2]; quadric4x4[0][3] = quadric4x4[3][0] = quadric[3]; quadric4x4[1][1] = quadric[4]; quadric4x4[1][2] = quadric4x4[2][1] = quadric[5]; quadric4x4[1][3] = quadric4x4[3][1] = quadric[6]; quadric4x4[2][2] = quadric[7]; quadric4x4[2][3] = quadric4x4[3][2] = quadric[8]; quadric4x4[3][3] = 1; // arbitrary value x = binId % this->NumberOfDivisions[0]; y = (binId / this->NumberOfDivisions[0]) % this->NumberOfDivisions[1]; z = binId / this->SliceSize; cellBounds[0] = this->Bounds[0] + x * this->XBinSize; cellBounds[1] = this->Bounds[0] + (x+1) * this->XBinSize; cellBounds[2] = this->Bounds[2] + y * this->YBinSize; cellBounds[3] = this->Bounds[2] + (y+1) * this->YBinSize; cellBounds[4] = this->Bounds[4] + z * this->ZBinSize; cellBounds[5] = this->Bounds[4] + (z+1) * this->ZBinSize; for (i = 0; i < 3; i++) { b[i] = -quadric4x4[3][i]; cellCenter[i] = (cellBounds[i*2+1] + cellBounds[i*2]) / 2.0; for (j = 0; j < 3; j++) { A[i][j] = quadric4x4[i][j]; } } // Compute diagonal matrix W // #define VTK_SVTHRESHOLD 1.0E-3 double maxW = 0.0, temp; vtkMath::SingularValueDecomposition3x3(A, U, w, VT); // Find maximum (magnitude) eigenvalue from SVD for (i = 0; i < 3; i++) { if ((temp = fabs(w[i])) > maxW) { maxW = temp; } } // Initialize matrix for (i = 0; i < 3; i++) { for (j = 0; j < 3; j++) { if (i == j) { if ( fabs(w[i] / maxW) > VTK_SVTHRESHOLD) { // If this is true, then w[i] != 0, so this division is ok. W[i][j] = 1.0/w[i]; } else { W[i][j] = 0.0; } } else { W[i][j] = 0.0; } } } #undef VTK_SVTHRESHOLD vtkMath::Transpose3x3(U, UT); vtkMath::Transpose3x3(VT, V); vtkMath::Multiply3x3(W, UT, tempMatrix); vtkMath::Multiply3x3(V, tempMatrix, tempMatrix); vtkMath::Multiply3x3(A, cellCenter, tempVector); for (i = 0; i < 3; i++) { tempVector[i] = b[i] - tempVector[i]; } vtkMath::Multiply3x3(tempMatrix, tempVector, tempVector); // Make absolutely sure that the point lies in the vicinity (enclosing // sphere) of the bin. If not, then clamp the point to the enclosing sphere. double deltaMag = vtkMath::Norm(tempVector); double radius = sqrt(this->XBinSize*this->XBinSize + this->YBinSize*this->YBinSize + this->ZBinSize*this->ZBinSize) / 2.0; if (deltaMag > radius) { tempVector[0] *= radius / deltaMag; tempVector[1] *= radius / deltaMag; tempVector[2] *= radius / deltaMag; } point[0] = cellCenter[0] + tempVector[0]; point[1] = cellCenter[1] + tempVector[1]; point[2] = cellCenter[2] + tempVector[2]; } //---------------------------------------------------------------------------- void vtkQuadricClustering::SetNumberOfDivisions(int div0, int div1, int div2) { this->SetNumberOfXDivisions(div0); this->SetNumberOfYDivisions(div1); this->SetNumberOfZDivisions(div2); } //---------------------------------------------------------------------------- void vtkQuadricClustering::SetNumberOfXDivisions(int num) { if (this->NumberOfXDivisions == num && this->ComputeNumberOfDivisions == 0) { return; } if (num < 1) { vtkErrorMacro("You cannot use less than one division."); return; } this->Modified(); this->NumberOfXDivisions = num; this->ComputeNumberOfDivisions = 0; } //---------------------------------------------------------------------------- void vtkQuadricClustering::SetNumberOfYDivisions(int num) { if (this->NumberOfYDivisions == num && this->ComputeNumberOfDivisions == 0) { return; } if (num < 1) { vtkErrorMacro("You cannot use less than one division."); return; } this->Modified(); this->NumberOfYDivisions = num; this->ComputeNumberOfDivisions = 0; } //---------------------------------------------------------------------------- void vtkQuadricClustering::SetNumberOfZDivisions(int num) { if (this->NumberOfZDivisions == num && this->ComputeNumberOfDivisions == 0) { return; } if (num < 1) { vtkErrorMacro("You cannot use less than one division."); return; } this->Modified(); this->NumberOfZDivisions = num; this->ComputeNumberOfDivisions = 0; } //---------------------------------------------------------------------------- int *vtkQuadricClustering::GetNumberOfDivisions() { static int divs[3]; this->GetNumberOfDivisions(divs); return divs; } //---------------------------------------------------------------------------- void vtkQuadricClustering::GetNumberOfDivisions(int divs[3]) { divs[0] = this->NumberOfXDivisions; divs[1] = this->NumberOfYDivisions; divs[2] = this->NumberOfZDivisions; } //---------------------------------------------------------------------------- void vtkQuadricClustering::SetDivisionOrigin(double x, double y, double z) { if (this->ComputeNumberOfDivisions && this->DivisionOrigin[0] == x && this->DivisionOrigin[1] == y && this->DivisionOrigin[2] == z) { return; } this->Modified(); this->DivisionOrigin[0] = x; this->DivisionOrigin[1] = y; this->DivisionOrigin[2] = z; this->ComputeNumberOfDivisions = 1; } //---------------------------------------------------------------------------- void vtkQuadricClustering::SetDivisionSpacing(double x, double y, double z) { if (this->ComputeNumberOfDivisions && this->DivisionSpacing[0] == x && this->DivisionSpacing[1] == y && this->DivisionSpacing[2] == z) { return; } if ( x <= 0 ) { vtkErrorMacro( << "Spacing (x) should be larger than 0.0, setting to 1.0" ); x = 1.0; } if ( y <= 0 ) { vtkErrorMacro( << "Spacing (y) should be larger than 0.0, setting to 1.0" ); y = 1.0; } if ( z <= 0 ) { vtkErrorMacro( << "Spacing (z) should be larger than 0.0, setting to 1.0" ); z = 1.0; } this->Modified(); this->DivisionSpacing[0] = x; this->DivisionSpacing[1] = y; this->DivisionSpacing[2] = z; this->ComputeNumberOfDivisions = 1; } //---------------------------------------------------------------------------- void vtkQuadricClustering::EndAppendUsingPoints(vtkPolyData *input, vtkPolyData *output) { vtkIdType i; vtkIdType outPtId; vtkPoints *inputPoints; vtkPoints *outputPoints; vtkIdType numPoints, numBins; vtkIdType binId; double *minError, e, pt[3]; double *q; inputPoints = input->GetPoints(); if (inputPoints == NULL) { return; } // Check for misuse of the Append methods. if (this->OutputTriangleArray == NULL || this->OutputLines == NULL) { vtkDebugMacro("Missing Array: Did you call StartAppend?"); return; } // Clean up if ( this->PreventDuplicateCells ) { delete this->CellSet; this->CellSet = NULL; } outputPoints = vtkPoints::New(); // Prepare to copy point data to output output->GetPointData()-> CopyAllocate(input->GetPointData(), this->NumberOfBinsUsed); // Allocate and initialize an array to hold errors for each bin. numBins = this->NumberOfDivisions[0] * this->NumberOfDivisions[1] * this->NumberOfDivisions[2]; minError = new double[numBins]; for (i = 0; i < numBins; ++i) { minError[i] = VTK_DOUBLE_MAX; } // Loop through the input points. numPoints = inputPoints->GetNumberOfPoints(); for (i = 0; i < numPoints; ++i) { inputPoints->GetPoint(i, pt); binId = this->HashPoint(pt); outPtId = this->QuadricArray[binId].VertexId; // Sanity check. if (outPtId == -1) { // This condition happens when there are points in the input that are // not used in any triangles, and therefore are never added to the // 3D hash structure. continue; } // Compute the error for this point. Note: the constant term is ignored. // It will be the same for every point in this bin, and it // is not stored in the quadric array anyway. q = this->QuadricArray[binId].Quadric; e = q[0]*pt[0]*pt[0] + 2.0*q[1]*pt[0]*pt[1] + 2.0*q[2]*pt[0]*pt[2] + 2.0*q[3]*pt[0] + q[4]*pt[1]*pt[1] + 2.0*q[5]*pt[1]*pt[2] + 2.0*q[6]*pt[1] + q[7]*pt[2]*pt[2] + 2.0*q[8]*pt[2]; if (e < minError[binId]) { minError[binId] = e; outputPoints->InsertPoint(outPtId, pt); // Since this is the same point as the input point, copy point data here too. output->GetPointData()->CopyData(input->GetPointData(),i,outPtId); } } output->SetPolys(this->OutputTriangleArray); output->SetPoints(outputPoints); outputPoints->Delete(); this->OutputTriangleArray->Delete(); this->OutputTriangleArray = NULL; if (this->OutputLines->GetNumberOfCells() > 0) { output->SetLines(this->OutputLines); } this->OutputLines->Delete(); this->OutputLines = NULL; this->EndAppendVertexGeometry(input, output); if (this->QuadricArray) { delete [] this->QuadricArray; this->QuadricArray = NULL; } delete [] minError; } //---------------------------------------------------------------------------- // This is not a perfect implementation, because it does not determine // which vertex cell is the best for a bin. The first detected is used. void vtkQuadricClustering::EndAppendVertexGeometry(vtkPolyData *input, vtkPolyData *output) { vtkCellArray *inVerts, *outVerts; vtkIdType *tmp = NULL; int tmpLength = 0; int tmpIdx; double pt[3]; int j; vtkIdType *ptIds = 0; vtkIdType numPts = 0; vtkIdType outPtId; vtkIdType binId, cellId, outCellId; inVerts = input->GetVerts(); outVerts = vtkCellArray::New(); for (cellId=0, inVerts->InitTraversal(); inVerts->GetNextCell(numPts, ptIds); cellId++) { if (tmpLength < numPts) { if (tmp) { delete[] tmp; } tmp = new vtkIdType[numPts]; tmpLength = numPts; } tmpIdx = 0; for (j = 0; j < numPts; ++j) { input->GetPoint(ptIds[j], pt); binId = this->HashPoint(pt); outPtId = this->QuadricArray[binId].VertexId; if (outPtId >= 0) { // Do not use this point. Destroy infomration in Quadric array. this->QuadricArray[binId].VertexId = -1; tmp[tmpIdx] = outPtId; ++tmpIdx; } } if (tmpIdx > 0) { // add poly vertex to output. outCellId = outVerts->InsertNextCell(tmpIdx, tmp); output->GetCellData()-> CopyData(input->GetCellData(), cellId, outCellId); } } if (tmp) { delete [] tmp; } if (outVerts->GetNumberOfCells() > 0) { output->SetVerts(outVerts); } outVerts->Delete(); } //---------------------------------------------------------------------------- // This method is called after the execution, but before the vertex array // is deleted. It changes some points to be based on the boundary edges. void vtkQuadricClustering::AppendFeatureQuadrics(vtkPolyData *pd, vtkPolyData *output) { vtkPolyData *input = vtkPolyData::New(); vtkPoints *edgePts; vtkCellArray *edges; vtkIdType i; vtkIdType binId; double featurePt[3]; // Find the boundary edges. input->ShallowCopy(pd); this->FeatureEdges->SetInput(input); this->FeatureEdges->Update(); edgePts = this->FeatureEdges->GetOutput()->GetPoints(); edges = this->FeatureEdges->GetOutput()->GetLines(); if (edges && edges->GetNumberOfCells() && edgePts) { this->AddEdges(edges, edgePts, 0, pd, output); if (this->UseFeaturePoints) { this->FindFeaturePoints(edges, edgePts, this->FeaturePointsAngle); for (i = 0; i < this->FeaturePoints->GetNumberOfPoints(); i++) { this->FeaturePoints->GetPoint(i, featurePt); binId = this->HashPoint(featurePt); this->AddVertex(binId, featurePt, 0, input, output); } } } // Release data. this->FeatureEdges->SetInputConnection(0, 0); this->FeatureEdges->GetOutput()->ReleaseData(); input->Delete(); } // Find the feature points of a given set of edges. // The points returned are (1) those used by only one edge, (2) those // used by > 2 edges, and (3) those where the angle between 2 edges // using this point is < angle. void vtkQuadricClustering::FindFeaturePoints(vtkCellArray *edges, vtkPoints *edgePts, double vtkNotUsed(angle)) { vtkIdType i, pointIds[2]; int j; vtkIdType *cellPts = 0; vtkIdType numCellPts; vtkIdList *pointIdList = vtkIdList::New(); vtkIdType numPts = edgePts->GetNumberOfPoints(); vtkIdType numEdges = edges->GetNumberOfCells(); vtkIdType edgeCount; vtkIdType **pointTable = new vtkIdType *[numPts]; double featurePoint[3]; double featureEdges[2][3]; double point1[3], point2[3]; vtkIdType *cellPointIds; double radAngle = vtkMath::RadiansFromDegrees( this->FeaturePointsAngle ); this->FeaturePoints->Allocate(numPts); for (i = 0; i < numPts; i++) { pointTable[i] = new vtkIdType[4]; pointTable[i][1] = 0; } edges->InitTraversal(); for (i = 0; i < numEdges; i++) { edges->GetNextCell(numCellPts, cellPts); for (j = 0; j < 2; j++) { pointIds[j] = pointIdList->InsertUniqueId(cellPts[j]); pointTable[pointIds[j]][0] = cellPts[j]; edgeCount = pointTable[pointIds[j]][1]; if (edgeCount < 2) { pointTable[pointIds[j]][edgeCount+2] = i; } pointTable[pointIds[j]][1]++; } } for (i = 0; i < numPts; i++) { if (pointTable[i][1] == 1) { edgePts->GetPoint(pointTable[i][0], featurePoint); } else if (pointTable[i][1] > 2) { edgePts->GetPoint(pointTable[i][0], featurePoint); } else if (pointTable[i][1] == 2) { for (j = 0; j < 2; j++) { edges->GetCell(3*pointTable[i][j+2], numCellPts, cellPointIds); if (cellPointIds[0] == pointTable[i][0]) { edgePts->GetPoint(cellPointIds[0], point1); edgePts->GetPoint(cellPointIds[1], point2); } else { edgePts->GetPoint(cellPointIds[1], point1); edgePts->GetPoint(cellPointIds[0], point2); } featureEdges[j][0] = point2[0] - point1[0]; featureEdges[j][1] = point2[1] - point1[1]; featureEdges[j][2] = point2[2] - point1[2]; vtkMath::Normalize(featureEdges[j]); } if (acos(vtkMath::Dot(featureEdges[0], featureEdges[1])) < radAngle) { edgePts->GetPoint(pointTable[i][0], featurePoint); } } } pointIdList->Delete(); for (i = 0; i < numPts; i++) { delete [] pointTable[i]; } delete [] pointTable; } //---------------------------------------------------------------------------- int vtkQuadricClustering::FillInputPortInformation(int port, vtkInformation *info) { int retval = this->Superclass::FillInputPortInformation(port, info); info->Set(vtkAlgorithm::INPUT_IS_OPTIONAL(), 1); return retval; } //---------------------------------------------------------------------------- void vtkQuadricClustering::PrintSelf(ostream& os, vtkIndent indent) { this->Superclass::PrintSelf(os,indent); os << indent << "Bounds: " << this->Bounds[0] << " " << this->Bounds[1] << " " << this->Bounds[2] << " " << this->Bounds[3] << " " << this->Bounds[4] << " " << this->Bounds[5] << "\n"; os << indent << "Use Input Points: " << (this->UseInputPoints ? "On\n" : "Off\n"); if (this->ComputeNumberOfDivisions) { os << indent << "Using Spacing and Origin to construct bins\n"; } else { os << indent << "Using input bounds and NumberOfDivisions to construct bins\n"; } os << indent << "Division Spacing: " << this->DivisionSpacing[0] << ", " << this->DivisionSpacing[1] << ", " << this->DivisionSpacing[2] << endl; os << indent << "Division Origin: " << this->DivisionOrigin[0] << ", " << this->DivisionOrigin[1] << ", " << this->DivisionOrigin[2] << endl; os << indent << "Number of X Divisions: " << this->NumberOfXDivisions << "\n"; os << indent << "Number of Y Divisions: " << this->NumberOfYDivisions << "\n"; os << indent << "Number of Z Divisions: " << this->NumberOfZDivisions << "\n"; os << indent << "Auto Adjust Number Of Divisions: " << (this->AutoAdjustNumberOfDivisions ? "On\n" : "Off\n"); os << indent << "Use Internal Triangles: " << (this->UseInternalTriangles ? "On\n" : "Off\n"); os << indent << "Use Feature Edges: " << this->UseFeatureEdges << endl; os << indent << "FeatureEdges: (" << this->FeatureEdges << ")\n"; os << indent << "Feature Points Angle: " << this->FeaturePointsAngle << endl; os << indent << "Use Feature Points: " << (this->UseFeaturePoints ? "On\n" : "Off\n"); os << indent << "Copy Cell Data : " << this->CopyCellData << endl; os << indent << "Prevent Duplicate Cells : " << (this->PreventDuplicateCells ? "On\n" : "Off\n"); }