diff --git a/vtkParallelTransportFrame.cxx b/vtkParallelTransportFrame.cxx index d459c76..c0f22df 100644 --- a/vtkParallelTransportFrame.cxx +++ b/vtkParallelTransportFrame.cxx @@ -36,6 +36,7 @@ vtkParallelTransportFrame::vtkParallelTransportFrame() this->SetTangentsArrayName("Tangents"); this->SetNormalsArrayName("Normals"); this->SetBinormalsArrayName("Binormals"); + this->SetRotationMinimizingFrames(false); } //---------------------------------------------------------------------------- @@ -119,7 +120,10 @@ int vtkParallelTransportFrame::RequestData( int numberOfCells = input->GetNumberOfCells(); for (int cellIndex = 0; cellIndex < numberOfCells; cellIndex++) { - this->ComputeAxisDirections(input, cellIndex, tangentsArray, normalsArray, binormalsArray); + if (this->RotationMinimizingFrames) + this->ComputeAxisDirections2(input, cellIndex, tangentsArray, normalsArray, binormalsArray); + else + this->ComputeAxisDirections(input, cellIndex, tangentsArray, normalsArray, binormalsArray); } output->GetPointData()->AddArray(tangentsArray); @@ -250,6 +254,165 @@ void vtkParallelTransportFrame::ComputeAxisDirections(vtkPolyData* input, vtkIdT } } +//---------------------------------------------------------------------------- +void vtkParallelTransportFrame::ComputeAxisDirections2(vtkPolyData* input, vtkIdType cellIndex, vtkDoubleArray* tangentsArray, vtkDoubleArray* normalsArray, vtkDoubleArray* binormalsArray) +{ + vtkPolyLine* polyLine = vtkPolyLine::SafeDownCast(input->GetCell(cellIndex)); + if (!polyLine) + { + return; + } + vtkIdType numberOfPointsInCell = polyLine->GetNumberOfPoints(); + if (numberOfPointsInCell < 2) + { + return; + } + + double tangent0[3] = { 0.0, 0.0, 0.0 }; + vtkIdType pointId0 = polyLine->GetPointId(0); + double pointPosition0[3]; + input->GetPoint(pointId0, pointPosition0); + + // Find tangent by direction vector by moving a minimal distance from the initial point + for (int pointIndex = 1; pointIndex < numberOfPointsInCell; pointIndex++) + { + vtkIdType pointId1 = polyLine->GetPointId(pointIndex); + double pointPosition1[3]; + input->GetPoint(pointId1, pointPosition1); + tangent0[0] = pointPosition1[0] - pointPosition0[0]; + tangent0[1] = pointPosition1[1] - pointPosition0[1]; + tangent0[2] = pointPosition1[2] - pointPosition0[2]; + if (vtkMath::Norm(tangent0) >= this->MinimumDistance) + { + break; + } + } + vtkMath::Normalize(tangent0); + + // Compute initial normal and binormal directions from the initial tangent and preferred + // normal/binormal directions. + double normal0[3] = {0.0, 0.0, 0.0}; + double binormal0[3] = {0.0, 0.0, 0.0}; + vtkMath::Cross(tangent0, this->PreferredInitialNormalVector, binormal0); + if (vtkMath::Norm(binormal0) > this->Tolerance) + { + vtkMath::Normalize(binormal0); + vtkMath::Cross(binormal0, tangent0, normal0); + } + else + { + vtkMath::Cross(this->PreferredInitialBinormalVector, tangent0, normal0); + vtkMath::Normalize(normal0); + vtkMath::Cross(tangent0, normal0, binormal0); + } + + tangentsArray->SetTuple(pointId0, tangent0); + normalsArray->SetTuple(pointId0, normal0); + binormalsArray->SetTuple(pointId0, binormal0); + + + + + + vtkIdType pointId2 = -1; + double tangent1[3] = { tangent0[0], tangent0[1], tangent0[2] }; + for (int i = 1; i < numberOfPointsInCell - 1; i++) + { + vtkIdType pointId1 = polyLine->GetPointId(i); + pointId2 = polyLine->GetPointId(i+1); + double pointPosition1[3]; + double pointPosition2[3]; + input->GetPoint(pointId1, pointPosition1); + input->GetPoint(pointId2, pointPosition2); + + tangent1[0] = pointPosition2[0] - pointPosition1[0]; + tangent1[1] = pointPosition2[1] - pointPosition1[1]; + tangent1[2] = pointPosition2[2] - pointPosition1[2]; + + vtkMath::Normalize(tangent1); + + tangentsArray->SetTuple(pointId1, tangent1); + + // Save current data for next iteration + tangent0[0] = tangent1[0]; + tangent0[1] = tangent1[1]; + tangent0[2] = tangent1[2]; + } + + if (pointId2 >= 0) + tangentsArray->SetTuple(pointId2, tangent1); + + + + vtkIdType pointIdi = -1; + vtkIdType pointIdip1 = -1; + double ri[3] = {normal0[0],normal0[1],normal0[2]}; + for (int i = 0; i < numberOfPointsInCell - 1; i++) + { + pointIdi = polyLine->GetPointId(i); + pointIdip1 = polyLine->GetPointId(i+1); + // Get positions + double xi[3]; + double xip1[3]; + input->GetPoint(pointIdi, xi); + input->GetPoint(pointIdip1, xip1); + // Get tangent vectors + double ti[3]; + double tip1[3]; + tangentsArray->GetTuple(pointIdi, ti); + tangentsArray->GetTuple(pointIdip1, tip1); + + // Compute reflection vector of R1 + double v1[3], v1_[3] = {0,0,0}; + v1[0] = v1_[0] = xip1[0] - xi[0]; + v1[1] = v1_[1] = xip1[1] - xi[1]; + v1[2] = v1_[2] = xip1[2] - xi[2]; + double c1 = 0; + c1 = vtkMath::Dot(v1,v1); + + // Compute rLi = R1*ri + double v1ri = 0; + v1ri = vtkMath::Dot(v1,ri); + vtkMath::MultiplyScalar(v1,v1ri); + vtkMath::MultiplyScalar(v1,-2/c1); + double rLi[3] = {0,0,0}; + vtkMath::Add(ri,v1,rLi); + + // Compute tLi = R1*ti + double v1ti = 0; + v1ti = vtkMath::Dot(v1_,ti); + vtkMath::MultiplyScalar(v1_,v1ti); + vtkMath::MultiplyScalar(v1_,-2/c1); + double tLi[3] = {0,0,0}; + vtkMath::Add(ti,v1_,tLi); + + // Compute reflection vector of R2 + double v2[3] = {0,0,0}; + vtkMath::Subtract(tip1,tLi,v2); + + double c2 = 0; + c2 = vtkMath::Dot(v2,v2); + + // Compute rip1 = R2*rLi (normal vector) + double v2rLi = 0; + v2rLi = vtkMath::Dot(v2,rLi); + vtkMath::MultiplyScalar(v2,v2rLi); + vtkMath::MultiplyScalar(v2,-2/c2); + double rip1[3] = {0,0,0}; + vtkMath::Add(rLi,v2,rip1); + vtkMath::Normalize(rip1); + + // Compute vector sip1 of the frame i+1 (binormal vector) + double sip1[3] = {0,0,0}; + vtkMath::Cross(tip1,rip1,sip1); + vtkMath::Normalize(sip1); + + //tangentsArray->SetTuple(pointIdip1, tip1); + normalsArray->SetTuple(pointIdip1, rip1); + binormalsArray->SetTuple(pointIdip1, sip1); + } +} + //---------------------------------------------------------------------------- void vtkParallelTransportFrame::RotateVector(double* inVector, double* outVector, const double* axis, double angle) { diff --git a/vtkParallelTransportFrame.h b/vtkParallelTransportFrame.h index c64cf78..6311866 100644 --- a/vtkParallelTransportFrame.h +++ b/vtkParallelTransportFrame.h @@ -26,10 +26,12 @@ /// /// References: /// - Parallel transport theory: R. Bishop, "There is more than one way to frame a curve", -/// American Mathematical Monthly, vol. 82, no. 3, pp. 246–251, 1975 +/// American Mathematical Monthly, vol. 82, no. 3, pp. 246�251, 1975 /// - Parallel transport implementation: Piccinelli M, Veneziani A, Steinman DA, Remuzzi A, Antiga L. /// "A framework for geometric analysis of vascular structures: application to cerebral aneurysms.", /// IEEE Trans Med Imaging. 2009 Aug;28(8):1141-55. doi: 10.1109/TMI.2009.2021652. +/// - Wang, W., J ̈uttler, B., Zheng, D., and Liu, Y. 2008. Computation of rotation minimizing frame. ACM Trans. Graph. 27, 1, Article 2 (March 2008), 18 pages. +/// DOI = 10.1145/1330511.1330513 http://doi.acm.org/10.1145/1330511.1330513 /// /// The initial implementation was based on VMTK (vtkvmtkCenterlineAttributesFilter) which was optimized /// and enhanced with more predictable initial normal vector direction. In the future, support for closed @@ -70,6 +72,13 @@ class VTK_ADDON_EXPORT vtkParallelTransportFrame : public vtkPolyDataAlgorithm vtkGetStringMacro(BinormalsArrayName); ///@} + ///@{ + /// Use rotation minimizing frames, otherwise use Bishop frames + /// By default is False + vtkSetMacro(RotationMinimizingFrames, vtkTypeBool); + vtkGetMacro(RotationMinimizingFrames, vtkTypeBool); + ///@} + /// Define the preferred direction of the normal vector at the first point of the curve. /// It is just "preferred" because the direction has to be orhogonal to the tangent, /// so in general the normal vector cannot point into exactly to a required direction. @@ -95,6 +104,7 @@ class VTK_ADDON_EXPORT vtkParallelTransportFrame : public vtkPolyDataAlgorithm int RequestData(vtkInformation *, vtkInformationVector **, vtkInformationVector *) override; void ComputeAxisDirections(vtkPolyData* input, vtkIdType cellIndex, vtkDoubleArray* tangentsArray, vtkDoubleArray* normalsArray, vtkDoubleArray* binormalsArray); + void ComputeAxisDirections2(vtkPolyData* input, vtkIdType cellIndex, vtkDoubleArray* tangentsArray, vtkDoubleArray* normalsArray, vtkDoubleArray* binormalsArray); /// Rotate a vector around an axis static void RotateVector(double* inVector, double* outVector, const double* axis, double angle); @@ -107,6 +117,8 @@ class VTK_ADDON_EXPORT vtkParallelTransportFrame : public vtkPolyDataAlgorithm char* NormalsArrayName = nullptr; char* BinormalsArrayName = nullptr; + vtkTypeBool RotationMinimizingFrames = false; + /// Tolerance value used for checking that a value is non-zero. double Tolerance = 1e-6; /// Minimum distance for comuting initial tangent direction