VTK实现局部平滑LocalSmoothing

/*========================================================================= Program: Visualization Toolkit Module: zxLocalSmoothingFilter.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 "zxLocalSmoothingFilter.h" #include "vtkCellArray.h" #include "vtkCellData.h" #include "vtkCellLocator.h" #include "vtkFloatArray.h" #include "vtkInformation.h" #include "vtkInformationVector.h" #include "vtkMath.h" #include "vtkObjectFactory.h" #include "vtkPointData.h" #include "vtkPolyData.h" #include "vtkPolygon.h" #include "vtkStreamingDemandDrivenPipeline.h" #include "vtkTriangleFilter.h" #include <limits> vtkStandardNewMacro(zxLocalSmoothingFilter); // The following code defines a helper class for performing mesh smoothing // across the surface of another mesh. typedef struct _vtkSmoothPoint { vtkIdType cellId; // cell int subId; // cell sub id double p[3]; // parametric coords in cell } vtkSmoothPoint; class vtkSmoothPoints { //;prevent man page generation public: vtkSmoothPoints(); ~vtkSmoothPoints() { delete[] this->Array; } vtkIdType GetNumberOfPoints() { return this->MaxId + 1; } vtkSmoothPoint* GetSmoothPoint(vtkIdType i) { return this->Array + i; } vtkSmoothPoint* InsertSmoothPoint(vtkIdType ptId) { if (ptId >= this->Size) { this->Resize(ptId + 1); } if (ptId > this->MaxId) { this->MaxId = ptId; } return this->Array + ptId; } vtkSmoothPoint* Resize(vtkIdType sz); // reallocates data void Reset() { this->MaxId = -1; } vtkSmoothPoint* Array; // pointer to data vtkIdType MaxId; // maximum index inserted thus far vtkIdType Size; // allocated size of data vtkIdType Extend; // grow array by this amount }; vtkSmoothPoints::vtkSmoothPoints() { this->MaxId = -1; this->Array = new vtkSmoothPoint[1000]; this->Size = 1000; this->Extend = 5000; } vtkSmoothPoint* vtkSmoothPoints::Resize(vtkIdType sz) { vtkSmoothPoint* newArray; vtkIdType newSize; if (sz >= this->Size) { newSize = this->Size + this->Extend * (((sz - this->Size) / this->Extend) + 1); } else { newSize = sz; } newArray = new vtkSmoothPoint[newSize]; memcpy(newArray, this->Array, (sz < this->Size ? sz : this->Size) * sizeof(vtkSmoothPoint)); this->Size = newSize; delete[] this->Array; this->Array = newArray; return this->Array; } // The following code defines methods for the zxLocalSmoothingFilter class // // Construct object with number of iterations 20; relaxation factor .01; // feature edge smoothing turned off; feature // angle 45 degrees; edge angle 15 degrees; and boundary smoothing turned // on. Error scalars and vectors are not generated (by default). The // convergence criterion is 0.0 of the bounding box diagonal. zxLocalSmoothingFilter::zxLocalSmoothingFilter() { this->Convergence = 0.0; // goes to number of specified iterations this->NumberOfIterations = 20; this->RelaxationFactor = .01; this->FeatureAngle = 45.0; this->EdgeAngle = 15.0; this->FeatureEdgeSmoothing = 0; this->BoundarySmoothing = 1; this->GenerateErrorScalars = 0; this->GenerateErrorVectors = 0; this->OutputPointsPrecision = vtkAlgorithm::DEFAULT_PRECISION; this->SmoothPoints = nullptr; this->m_smoothingIds = vtkIdList::New(); // optional second input this->SetNumberOfInputPorts(2); } void zxLocalSmoothingFilter::SetSourceData(vtkPolyData* source) { this->SetInputData(1, source); } vtkPolyData* zxLocalSmoothingFilter::GetSource() { if (this->GetNumberOfInputConnections(1) < 1) { return nullptr; } return vtkPolyData::SafeDownCast(this->GetExecutive()->GetInputData(1, 0)); } #define VTK_SIMPLE_VERTEX 0 #define VTK_FIXED_VERTEX 1 #define VTK_FEATURE_EDGE_VERTEX 2 #define VTK_BOUNDARY_EDGE_VERTEX 3 namespace { // Special structure for marking vertices typedef struct _vtkMeshVertex { char type; vtkIdList* edges; // connected edges (list of connected point ids) _vtkMeshVertex() { type = VTK_FIXED_VERTEX; // can smooth edges = nullptr; } } vtkMeshVertex, * vtkMeshVertexPtr; template <typename T> struct vtkSPDF_InternalParams { zxLocalSmoothingFilter* spdf; int numberOfIterations; vtkPoints* newPts; T factor; T conv; vtkIdType numPts; vtkMeshVertexPtr vertexPtr; vtkPolyData* source; vtkSmoothPoints* SmoothPoints; double* w; vtkCellLocator* cellLocator; }; template <typename T> void vtkSPDF_MovePoints(vtkSPDF_InternalParams<T>& params) { int iterationNumber = 0; for (T maxDist = std::numeric_limits<T>::max(); maxDist > params.conv && iterationNumber < params.numberOfIterations; ++iterationNumber) { if (iterationNumber && !(iterationNumber % 5)) { params.spdf->UpdateProgress(0.5 + 0.5 * iterationNumber / params.numberOfIterations); if (params.spdf->GetAbortExecute()) { break; } } maxDist = 0.0; T* newPtsCoords = static_cast<T*>(params.newPts->GetVoidPointer(0)); T* start = newPtsCoords; vtkMeshVertexPtr vertsPtr = params.vertexPtr; vtkIdType npts, * edgeIdPtr; T dist, deltaX[3]; double dist2, xNew[3], closestPt[3]; // For each non-fixed vertex of the mesh, move the point toward the mean // position of its connected neighbors using the relaxation factor. for (vtkIdType i = 0; i < params.numPts; ++i) { if (vertsPtr->type != VTK_FIXED_VERTEX && vertsPtr->edges != nullptr && (npts = vertsPtr->edges->GetNumberOfIds()) > 0) { deltaX[0] = deltaX[1] = deltaX[2] = 0.0; edgeIdPtr = vertsPtr->edges->GetPointer(0); // Compute the mean (cumulated) direction vector for (vtkIdType j = 0; j < npts; ++j) { for (unsigned short k = 0; k < 3; ++k) { deltaX[k] += *(start + 3 * (*edgeIdPtr) + k); } ++edgeIdPtr; } // for all connected points // Move the point *newPtsCoords += params.factor * (deltaX[0] / npts - (*newPtsCoords)); xNew[0] = *newPtsCoords; ++newPtsCoords; *newPtsCoords += params.factor * (deltaX[1] / npts - (*newPtsCoords)); xNew[1] = *newPtsCoords; ++newPtsCoords; *newPtsCoords += params.factor * (deltaX[2] / npts - (*newPtsCoords)); xNew[2] = *newPtsCoords; ++newPtsCoords; // Constrain point to surface if (params.source) { vtkSmoothPoint* sPtr = params.SmoothPoints->GetSmoothPoint(i