From 811b81a15501917135de75d5ece77b44312bdb57 Mon Sep 17 00:00:00 2001 From: "Mauro I. Dominguez" Date: Thu, 21 Jul 2022 18:28:44 -0300 Subject: [PATCH 1/3] ENH: add Rotation Minimizing Frames --- vtkParallelTransportFrame.cxx | 161 +++++++++++++++++++++++++++++++++- vtkParallelTransportFrame.h | 3 +- 2 files changed, 162 insertions(+), 2 deletions(-) diff --git a/vtkParallelTransportFrame.cxx b/vtkParallelTransportFrame.cxx index d459c76..1e8ec4f 100644 --- a/vtkParallelTransportFrame.cxx +++ b/vtkParallelTransportFrame.cxx @@ -119,7 +119,7 @@ int vtkParallelTransportFrame::RequestData( int numberOfCells = input->GetNumberOfCells(); for (int cellIndex = 0; cellIndex < numberOfCells; cellIndex++) { - this->ComputeAxisDirections(input, cellIndex, tangentsArray, normalsArray, binormalsArray); + this->ComputeAxisDirections2(input, cellIndex, tangentsArray, normalsArray, binormalsArray); } output->GetPointData()->AddArray(tangentsArray); @@ -250,6 +250,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..ce9aa50 100644 --- a/vtkParallelTransportFrame.h +++ b/vtkParallelTransportFrame.h @@ -26,7 +26,7 @@ /// /// 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. @@ -95,6 +95,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); From 401a1909cda70c58c2c1991d458fbd5e75cd4b4c Mon Sep 17 00:00:00 2001 From: "Mauro I. Dominguez" Date: Fri, 22 Jul 2022 06:51:19 -0300 Subject: [PATCH 2/3] ENH: add Rotation Minimizing Frames calculation --- vtkParallelTransportFrame.cxx | 6 +++++- vtkParallelTransportFrame.h | 9 +++++++++ 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/vtkParallelTransportFrame.cxx b/vtkParallelTransportFrame.cxx index 1e8ec4f..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->ComputeAxisDirections2(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); diff --git a/vtkParallelTransportFrame.h b/vtkParallelTransportFrame.h index ce9aa50..1d32699 100644 --- a/vtkParallelTransportFrame.h +++ b/vtkParallelTransportFrame.h @@ -70,6 +70,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. @@ -108,6 +115,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 From ad0feb9404ae217e25f84d51f2c338624b43d873 Mon Sep 17 00:00:00 2001 From: "Mauro I. Dominguez" Date: Fri, 22 Jul 2022 06:54:40 -0300 Subject: [PATCH 3/3] ENH: add reference to paper source of the RMF algorithm --- vtkParallelTransportFrame.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/vtkParallelTransportFrame.h b/vtkParallelTransportFrame.h index 1d32699..6311866 100644 --- a/vtkParallelTransportFrame.h +++ b/vtkParallelTransportFrame.h @@ -30,6 +30,8 @@ /// - 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