Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Rotation minimizing frames #28

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Rotation minimizing frames #28

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
811b81a
ENH: add Rotation Minimizing Frames
mauigna06 Jul 21, 2022
401a190
ENH: add Rotation Minimizing Frames calculation
mauigna06 Jul 22, 2022
ad0feb9
ENH: add reference to paper source of the RMF algorithm
mauigna06 Jul 22, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
165 changes: 164 additions & 1 deletion vtkParallelTransportFrame.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@ vtkParallelTransportFrame::vtkParallelTransportFrame()
this->SetTangentsArrayName("Tangents");
this->SetNormalsArrayName("Normals");
this->SetBinormalsArrayName("Binormals");
this->SetRotationMinimizingFrames(false);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We initialize simple types in the .h file. Please do it for this method selection variable, too. (it will be an enum not a bool, see below).

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Parallel transport, Double reflection method, etc. all minimize the rotation. The method name in the Wang2008 paper is "double reflection method". Please rename this class to vtkRotationMinimizingFrame and add a Method enum property (follow the VTK conventions for Get/Set method naming). Method names will be ParallelTransport and DoubleReflection. Add a deprecated vtkParallelTransfportFrame class, derived from vtkRotationMinimizingFrame (it does not add or remove any feature, just makes the vtkParallelTransfportFrame available for a while for backward compatibility).

}

//----------------------------------------------------------------------------
Expand Down Expand Up @@ -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);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

According to VTK style, braces {} are always required in branches, even if there is a single instruction in the branch. Please add them.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Rename the methods to ComputeAxisDirectionsParallelTransport and ComputeAxisDirectionsDoubleReflection.

else
this->ComputeAxisDirections(input, cellIndex, tangentsArray, normalsArray, binormalsArray);
}

output->GetPointData()->AddArray(tangentsArray);
Expand Down Expand Up @@ -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];
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Always initialize vectors to avoid having different behavior in debug and release mode (I see that this was the original implementation, so change it there, too).

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++)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is just too much copy-paste from ComputeAxisDirections. Factor our the common parts into separate methods and call them from both ComputeAxisDirections variants.

{
vtkIdType pointId1 = polyLine->GetPointId(pointIndex);
double pointPosition1[3];
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Initialize to {0.0, 0.0, 0.0} unless you find that there is a measurable performance penalty.

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);


Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove the extra newlines. Add a comment if you want to indicate a new section of the code.




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);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

missing braces



Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

remove extra newlines


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);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

uncomment or remove the line

normalsArray->SetTuple(pointIdip1, rip1);
binormalsArray->SetTuple(pointIdip1, sip1);
}
}

//----------------------------------------------------------------------------
void vtkParallelTransportFrame::RotateVector(double* inVector, double* outVector, const double* axis, double angle)
{
Expand Down
14 changes: 13 additions & 1 deletion vtkParallelTransportFrame.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We only use ASCII characters in the code. Please replace the unicode character by simple hyphen -.

/// - 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.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please describe the most important characteristics of the two methods (especially when one should be used vs. the other), do a simple performance profiling (in release mode, compute directions for many large curves and see if the difference is noticeable).

/// DOI = 10.1145/1330511.1330513 http://doi.acm.org/10.1145/1330511.1330513
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Add open-access URL as well (https://www.microsoft.com/en-us/research/wp-content/uploads/2016/12/Computation-of-rotation-minimizing-frames.pdf)

///
/// 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
Expand Down Expand Up @@ -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.
Expand All @@ -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);
Expand All @@ -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
Expand Down