+/*=========================================================================\r
+ *\r
+ * Copyright NumFOCUS\r
+ *\r
+ * Licensed under the Apache License, Version 2.0 (the "License");\r
+ * you may not use this file except in compliance with the License.\r
+ * You may obtain a copy of the License at\r
+ *\r
+ * http://www.apache.org/licenses/LICENSE-2.0.txt\r
+ *\r
+ * Unless required by applicable law or agreed to in writing, software\r
+ * distributed under the License is distributed on an "AS IS" BASIS,\r
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\r
+ * See the License for the specific language governing permissions and\r
+ * limitations under the License.\r
+ *\r
+ *=========================================================================*/\r
+\r
+#ifndef itkFastMarchingQuadEdgeMeshFilterResultsBase_hxx\r
+#define itkFastMarchingQuadEdgeMeshFilterResultsBase_hxx\r
+\r
+#include "itkFastMarchingQuadEdgeMeshFilterResultsBase.h"\r
+\r
+#include "itkMath.h"\r
+#include "itkVector.h"\r
+#include "itkMatrix.h"\r
+\r
+namespace itk\r
+{\r
+template <typename TInput, typename TOutput>\r
+FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::FastMarchingQuadEdgeMeshFilterResultsBase()\r
+ : Superclass()\r
+{\r
+ this->m_InputMesh = nullptr;\r
+ this->m_BackInputMesh = nullptr;\r
+}\r
+\r
+template <typename TInput, typename TOutput>\r
+IdentifierType\r
+FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::GetTotalNumberOfNodes() const\r
+{\r
+ return this->GetInput()->GetNumberOfPoints();\r
+}\r
+\r
+template <typename TInput, typename TOutput>\r
+void\r
+FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::SetOutputValue(OutputMeshType * oMesh,\r
+ const NodeType & iNode,\r
+ const OutputPixelType & iValue)\r
+{\r
+ oMesh->SetPointData(iNode, iValue);\r
+}\r
+\r
+template <typename TInput, typename TOutput>\r
+const typename FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::OutputPixelType\r
+FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::GetOutputValue(OutputMeshType * oMesh,\r
+ const NodeType & iNode) const\r
+{\r
+ OutputPixelType outputValue = NumericTraits<OutputPixelType>::ZeroValue();\r
+ oMesh->GetPointData(iNode, &outputValue);\r
+ return outputValue;\r
+}\r
+\r
+\r
+template <typename TInput, typename TOutput>\r
+unsigned char\r
+FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::GetLabelValueForGivenNode(const NodeType & iNode) const\r
+{\r
+ auto it = m_Label.find(iNode);\r
+\r
+ if (it != m_Label.end())\r
+ {\r
+ return it->second;\r
+ }\r
+ else\r
+ {\r
+ return Traits::Far;\r
+ }\r
+}\r
+\r
+template <typename TInput, typename TOutput>\r
+void\r
+FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::SetLabelValueForGivenNode(const NodeType & iNode,\r
+ const LabelType & iLabel)\r
+{\r
+ m_Label[iNode] = iLabel;\r
+}\r
+\r
+template <typename TInput, typename TOutput>\r
+void\r
+FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::UpdateNeighbors(OutputMeshType * oMesh, const NodeType & iNode)\r
+{\r
+ //printf("PG FastMarchingQuadEdgeMeshFilterResultsBase::UpdateNeighbors Start \n");\r
+ OutputPointType p;\r
+ oMesh->GetPoint(iNode, &p);\r
+\r
+ OutputQEType * qe = p.GetEdge();\r
+\r
+ if (qe)\r
+ {\r
+ OutputQEType * qe_it = qe;\r
+ this->m_InputMesh = this->GetInput();\r
+ do\r
+ {\r
+ if (qe_it)\r
+ {\r
+ OutputPointIdentifierType neigh_id = qe_it->GetDestination();\r
+\r
+ const char label = this->GetLabelValueForGivenNode(neigh_id);\r
+\r
+ if ((label != Traits::Alive) && (label != Traits::InitialTrial) && (label != Traits::Forbidden))\r
+ {\r
+ this->UpdateValue(oMesh, neigh_id);\r
+ }\r
+ }\r
+ else\r
+ {\r
+ itkGenericExceptionMacro(<< "qe_it is nullptr");\r
+ }\r
+ qe_it = qe_it->GetOnext();\r
+ } while (qe_it != qe);\r
+ }\r
+ else\r
+ {\r
+ itkGenericExceptionMacro(<< "qe is nullptr");\r
+ }\r
+ //printf("PG FastMarchingQuadEdgeMeshFilterResultsBase::UpdateNeighbors End \n");\r
+}\r
+\r
+template <typename TInput, typename TOutput>\r
+void\r
+FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::UpdateValue(OutputMeshType * oMesh, const NodeType & iNode)\r
+{\r
+ OutputPointType p;\r
+ oMesh->GetPoint(iNode, &p);\r
+\r
+ InputPixelType F = NumericTraits<InputPixelType>::ZeroValue();\r
+ //this->m_InputMesh->GetPointData(iNode, &F);\r
+ F = 1.;\r
+ if (F < 0.)\r
+ {\r
+ itkGenericExceptionMacro(<< "F < 0.");\r
+ }\r
+\r
+ OutputQEType * qe = p.GetEdge();\r
+\r
+ OutputPixelType outputPixel = this->m_LargeValue;\r
+\r
+ if (qe)\r
+ {\r
+ OutputQEType * qe_it = qe;\r
+ qe_it = qe_it->GetOnext();\r
+\r
+ do\r
+ {\r
+ OutputQEType * qe_it2 = qe_it->GetOnext();\r
+\r
+ if (qe_it2)\r
+ {\r
+ if (qe_it->GetLeft() != OutputMeshType::m_NoFace)\r
+ {\r
+ OutputPointIdentifierType id1 = qe_it->GetDestination();\r
+ OutputPointIdentifierType id2 = qe_it2->GetDestination();\r
+\r
+ const auto label1 = static_cast<LabelType>(this->GetLabelValueForGivenNode(id1));\r
+ const auto label2 = static_cast<LabelType>(this->GetLabelValueForGivenNode(id2));\r
+\r
+ bool IsFar1 = (label1 != Traits::Far);\r
+ bool IsFar2 = (label2 != Traits::Far);\r
+\r
+ if (IsFar1 || IsFar2)\r
+ {\r
+ OutputPointType q1 = oMesh->GetPoint(id1);\r
+ OutputPointType q2 = oMesh->GetPoint(id2);\r
+\r
+ auto val1 = static_cast<OutputVectorRealType>(this->GetOutputValue(oMesh, id1));\r
+ auto val2 = static_cast<OutputVectorRealType>(this->GetOutputValue(oMesh, id2));\r
+\r
+ if (val1 > val2)\r
+ {\r
+ OutputPointType temp_pt(q1);\r
+ q1 = q2;\r
+ q2 = temp_pt;\r
+\r
+ std::swap(val1, val2);\r
+ std::swap(IsFar1, IsFar2);\r
+ std::swap(id1, id2);\r
+ }\r
+\r
+ const OutputVectorRealType temp =\r
+ this->Solve(oMesh, iNode, p, F, id1, q1, IsFar1, val1, id2, q2, IsFar2, val2);\r
+\r
+ outputPixel = std::min(outputPixel, static_cast<OutputPixelType>(temp));\r
+ }\r
+ }\r
+\r
+ qe_it = qe_it2;\r
+ }\r
+ else\r
+ {\r
+ // throw one exception here\r
+ itkGenericExceptionMacro(<< "qe_it2 is nullptr");\r
+ }\r
+ } while (qe_it != qe);\r
+\r
+ if (outputPixel < this->m_LargeValue)\r
+ {\r
+ this->SetOutputValue(oMesh, iNode, outputPixel);\r
+\r
+ this->SetLabelValueForGivenNode(iNode, Traits::Trial);\r
+\r
+ this->m_Heap.push(NodePairType(iNode, outputPixel));\r
+ }\r
+ }\r
+ else\r
+ {\r
+ // throw one exception\r
+ itkGenericExceptionMacro(<< "qe_it is nullptr");\r
+ }\r
+}\r
+\r
+template <typename TInput, typename TOutput>\r
+const typename FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::OutputVectorRealType\r
+FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::Solve(OutputMeshType * oMesh,\r
+ const NodeType & iId,\r
+ const OutputPointType & iCurrentPoint,\r
+ const OutputVectorRealType & iF,\r
+ const NodeType & iId1,\r
+ const OutputPointType & iP1,\r
+ const bool & iIsFar1,\r
+ const OutputVectorRealType iVal1,\r
+ const NodeType & iId2,\r
+ const OutputPointType & iP2,\r
+ const bool & iIsFar2,\r
+ const OutputVectorRealType & iVal2) const\r
+{\r
+ OutputVectorType Edge1 = iP1 - iCurrentPoint;\r
+ OutputVectorType Edge2 = iP2 - iCurrentPoint;\r
+\r
+ OutputVectorRealType sq_norm1 = Edge1.GetSquaredNorm();\r
+ OutputVectorRealType norm1 = 0.;\r
+\r
+ OutputVectorRealType epsilon = NumericTraits<OutputVectorRealType>::epsilon();\r
+\r
+ if (sq_norm1 > epsilon)\r
+ {\r
+ norm1 = std::sqrt(sq_norm1);\r
+\r
+ OutputVectorRealType inv_norm1 = 1. / norm1;\r
+ Edge1 *= inv_norm1;\r
+ }\r
+\r
+ OutputVectorRealType sq_norm2 = Edge2.GetSquaredNorm();\r
+ OutputVectorRealType norm2 = 0.;\r
+\r
+ if (sq_norm2 > epsilon)\r
+ {\r
+ norm2 = std::sqrt(sq_norm2);\r
+\r
+ OutputVectorRealType inv_norm2 = 1. / norm2;\r
+ Edge2 *= inv_norm2;\r
+ }\r
+\r
+ if (!iIsFar1 && iIsFar2)\r
+ {\r
+ // only one point is a contributor\r
+ return iVal2 + norm2 * iF;\r
+ }\r
+ if (iIsFar1 && !iIsFar2)\r
+ {\r
+ // only one point is a contributor\r
+ return iVal1 + norm1 * iF;\r
+ }\r
+\r
+ if (iIsFar1 && iIsFar2)\r
+ {\r
+ auto dot = static_cast<OutputVectorRealType>(Edge1 * Edge2);\r
+\r
+ if (dot >= 0.)\r
+ {\r
+ return ComputeUpdate(iVal1, iVal2, norm2, sq_norm2, norm1, sq_norm1, dot, iF);\r
+ }\r
+ else\r
+ {\r
+ OutputVectorRealType sq_norm3, norm3, dot1, dot2;\r
+ OutputPointIdentifierType new_id;\r
+\r
+ bool unfolded =\r
+ UnfoldTriangle(oMesh, iId, iCurrentPoint, iId1, iP1, iId2, iP2, norm3, sq_norm3, dot1, dot2, new_id);\r
+\r
+ if (unfolded)\r
+ {\r
+ OutputVectorRealType t_sq_norm3 = norm3 * norm3;\r
+\r
+ auto val3 = static_cast<OutputVectorRealType>(this->GetOutputValue(oMesh, new_id));\r
+ OutputVectorRealType t1 = ComputeUpdate(iVal1, val3, norm3, t_sq_norm3, norm1, sq_norm1, dot1, iF);\r
+ OutputVectorRealType t2 = ComputeUpdate(iVal2, val3, norm3, t_sq_norm3, norm2, sq_norm2, dot2, iF);\r
+\r
+ return std::min(t1, t2);\r
+ }\r
+ else\r
+ {\r
+ return ComputeUpdate(iVal1, iVal2, norm2, sq_norm2, norm1, sq_norm1, dot, iF);\r
+ }\r
+ }\r
+ }\r
+\r
+ return this->m_LargeValue;\r
+}\r
+\r
+template <typename TInput, typename TOutput>\r
+const typename FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::OutputVectorRealType\r
+FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::ComputeUpdate(const OutputVectorRealType & iVal1,\r
+ const OutputVectorRealType & iVal2,\r
+ const OutputVectorRealType & iNorm1,\r
+ const OutputVectorRealType & iSqNorm1,\r
+ const OutputVectorRealType & iNorm2,\r
+ const OutputVectorRealType & iSqNorm2,\r
+ const OutputVectorRealType & iDot,\r
+ const OutputVectorRealType & iF) const\r
+{\r
+ auto large_value = static_cast<OutputVectorRealType>(this->m_LargeValue);\r
+\r
+ OutputVectorRealType t = large_value;\r
+\r
+ OutputVectorRealType CosAngle = iDot;\r
+ OutputVectorRealType SinAngle = std::sqrt(1. - iDot * iDot);\r
+\r
+ OutputVectorRealType u = iVal2 - iVal1;\r
+\r
+ OutputVectorRealType sq_u = u * u;\r
+ OutputVectorRealType f2 = iSqNorm1 + iSqNorm2 - 2. * iNorm1 * iNorm2 * CosAngle;\r
+ OutputVectorRealType f1 = iNorm2 * u * (iNorm1 * CosAngle - iNorm2);\r
+ OutputVectorRealType f0 = iSqNorm2 * (sq_u - iF * iF * iSqNorm1 * SinAngle * SinAngle);\r
+\r
+ OutputVectorRealType delta = f1 * f1 - f0 * f2;\r
+\r
+ OutputVectorRealType epsilon = NumericTraits<OutputVectorRealType>::epsilon();\r
+\r
+ if (delta >= 0.)\r
+ {\r
+ if (itk::Math::abs(f2) > epsilon)\r
+ {\r
+ t = (-f1 - std::sqrt(delta)) / f2;\r
+\r
+ // test if we must must choose the other solution\r
+ if ((t < u) || (iNorm2 * (t - u) / t < iNorm1 * CosAngle) || (iNorm1 / CosAngle < iNorm2 * (t - u) / t))\r
+ {\r
+ t = (-f1 + std::sqrt(delta)) / f2;\r
+ }\r
+ }\r
+ else\r
+ {\r
+ if (itk::Math::abs(f1) > epsilon)\r
+ {\r
+ t = -f0 / f1;\r
+ }\r
+ else\r
+ {\r
+ t = -large_value;\r
+ }\r
+ }\r
+ }\r
+ else\r
+ {\r
+ t = -large_value;\r
+ }\r
+\r
+ // choose the update from the 2 vertex only if upwind criterion is met\r
+ if ((u < t) && (iNorm1 * CosAngle < iNorm2 * (t - u) / t) && (iNorm2 * (t - u) / t < iNorm1 / CosAngle))\r
+ {\r
+ return t + iVal1;\r
+ }\r
+ else\r
+ {\r
+ return std::min(iNorm2 * iF + iVal1, iNorm1 * iF + iVal2);\r
+ }\r
+}\r
+\r
+template <typename TInput, typename TOutput>\r
+bool\r
+FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::UnfoldTriangle(OutputMeshType * oMesh,\r
+ const OutputPointIdentifierType & iId,\r
+ const OutputPointType & iP,\r
+ const OutputPointIdentifierType & iId1,\r
+ const OutputPointType & iP1,\r
+ const OutputPointIdentifierType & iId2,\r
+ const OutputPointType & iP2,\r
+ OutputVectorRealType & oNorm,\r
+ OutputVectorRealType & oSqNorm,\r
+ OutputVectorRealType & oDot1,\r
+ OutputVectorRealType & oDot2,\r
+ OutputPointIdentifierType & oId) const\r
+{\r
+ (void)iId;\r
+\r
+ OutputVectorType Edge1 = iP1 - iP;\r
+ OutputVectorRealType Norm1 = Edge1.GetNorm();\r
+\r
+ OutputVectorRealType epsilon = NumericTraits<OutputVectorRealType>::epsilon();\r
+\r
+ if (Norm1 > epsilon)\r
+ {\r
+ OutputVectorRealType inv_Norm = 1. / Norm1;\r
+ Edge1 *= inv_Norm;\r
+ }\r
+\r
+ OutputVectorType Edge2 = iP2 - iP;\r
+ OutputVectorRealType Norm2 = Edge2.GetNorm();\r
+ if (Norm2 > epsilon)\r
+ {\r
+ OutputVectorRealType inv_Norm = 1. / Norm2;\r
+ Edge2 *= inv_Norm;\r
+ }\r
+\r
+ auto dot = static_cast<OutputVectorRealType>(Edge1 * Edge2);\r
+\r
+ // the equation of the lines defining the unfolding region\r
+ // [e.g. line 1 : {x ; <x,eq1>=0} ]\r
+ using Vector2DType = Vector<OutputVectorRealType, 2>;\r
+ using Matrix2DType = Matrix<OutputVectorRealType, 2, 2>;\r
+\r
+ Vector2DType v1;\r
+ v1[0] = dot;\r
+ v1[1] = std::sqrt(1. - dot * dot);\r
+\r
+ Vector2DType v2;\r
+ v2[0] = 1.;\r
+ v2[1] = 0.;\r
+\r
+ Vector2DType x1;\r
+ x1[0] = Norm1;\r
+ x1[1] = 0.;\r
+\r
+ Vector2DType x2 = Norm2 * v1;\r
+\r
+ // keep track of the starting point\r
+ Vector2DType x_start1(x1);\r
+ Vector2DType x_start2(x2);\r
+\r
+ OutputPointIdentifierType id1 = iId1;\r
+ OutputPointIdentifierType id2 = iId2;\r
+\r
+ OutputQEType * qe = oMesh->FindEdge(id1, id2);\r
+ qe = qe->GetSym();\r
+\r
+ OutputPointType t_pt1 = iP1;\r
+ OutputPointType t_pt2 = iP2;\r
+ OutputPointType t_pt = t_pt1;\r
+\r
+ unsigned int nNum = 0;\r
+ while (nNum < 50 && qe->GetLeft() != OutputMeshType::m_NoFace)\r
+ {\r
+ OutputQEType * qe_Lnext = qe->GetLnext();\r
+ OutputPointIdentifierType t_id = qe_Lnext->GetDestination();\r
+\r
+ oMesh->GetPoint(t_id, &t_pt);\r
+\r
+ Edge1 = t_pt2 - t_pt1;\r
+ Norm1 = Edge1.GetNorm();\r
+ if (Norm1 > epsilon)\r
+ {\r
+ OutputVectorRealType inv_Norm = 1. / Norm1;\r
+ Edge1 *= inv_Norm;\r
+ }\r
+\r
+ Edge2 = t_pt - t_pt1;\r
+ Norm2 = Edge2.GetNorm();\r
+ if (Norm2 > epsilon)\r
+ {\r
+ OutputVectorRealType inv_Norm = 1. / Norm2;\r
+ Edge2 *= inv_Norm;\r
+ }\r
+\r
+ /* compute the position of the new point x on the unfolding plane (via a rotation of -alpha on (x2-x1)/rNorm1 )\r
+ | cos(alpha) sin(alpha)|\r
+ x = |-sin(alpha) cos(alpha)| * [x2-x1]*rNorm2/rNorm1 + x1 where cos(alpha)=dot\r
+ */\r
+ Vector2DType vv;\r
+ if (Norm1 > epsilon)\r
+ {\r
+ vv = (x2 - x1) * Norm2 / Norm1;\r
+ }\r
+ else\r
+ {\r
+ vv.Fill(0.);\r
+ }\r
+\r
+ dot = Edge1 * Edge2;\r
+\r
+ Matrix2DType rotation;\r
+ rotation[0][0] = dot;\r
+ rotation[0][1] = std::sqrt(1. - dot * dot);\r
+ rotation[1][0] = -rotation[0][1];\r
+ rotation[1][1] = dot;\r
+\r
+ Vector2DType x = rotation * vv + x1;\r
+\r
+\r
+ /* compute the intersection points.\r
+ We look for x=x1+lambda*(x-x1) or x=x2+lambda*(x-x2) with <x,eqi>=0, so */\r
+ OutputVectorRealType lambda11 = -(x1 * v1) / ((x - x1) * v1); // left most\r
+ OutputVectorRealType lambda12 = -(x1 * v2) / ((x - x1) * v2); // right most\r
+ OutputVectorRealType lambda21 = -(x2 * v1) / ((x - x2) * v1); // left most\r
+ OutputVectorRealType lambda22 = -(x2 * v2) / ((x - x2) * v2); // right most\r
+ bool bIntersect11 = (lambda11 >= 0.) && (lambda11 <= 1.);\r
+ bool bIntersect12 = (lambda12 >= 0.) && (lambda12 <= 1.);\r
+ bool bIntersect21 = (lambda21 >= 0.) && (lambda21 <= 1.);\r
+ bool bIntersect22 = (lambda22 >= 0.) && (lambda22 <= 1.);\r
+\r
+ if (bIntersect11 && bIntersect12)\r
+ {\r
+ qe = oMesh->FindEdge(id1, t_id);\r
+ qe = qe->GetSym();\r
+ id2 = t_id;\r
+ t_pt2 = t_pt;\r
+ x2 = x;\r
+ }\r
+ else\r
+ {\r
+ if (bIntersect21 && bIntersect22)\r
+ {\r
+ qe = oMesh->FindEdge(id2, t_id);\r
+ qe = qe->GetSym();\r
+ id1 = t_id;\r
+ t_pt1 = t_pt;\r
+ x1 = x;\r
+ }\r
+ else\r
+ { // bIntersect11 && !bIntersect12 && !bIntersect21 && bIntersect22\r
+\r
+ // that's it, we have found the point\r
+ oSqNorm = x.GetSquaredNorm();\r
+\r
+ if (oSqNorm > epsilon)\r
+ {\r
+ oNorm = std::sqrt(oSqNorm);\r
+ OutputVectorRealType temp_norm = x_start1.GetNorm();\r
+ if (temp_norm > epsilon)\r
+ {\r
+ oDot1 = x * x_start1 / (oNorm * temp_norm);\r
+ }\r
+ else\r
+ {\r
+ oDot1 = 0.;\r
+ }\r
+ temp_norm = x_start2.GetNorm();\r
+ if (temp_norm > epsilon)\r
+ {\r
+ oDot2 = x * x_start2 / (oNorm * temp_norm);\r
+ }\r
+ else\r
+ {\r
+ oDot2 = 0.;\r
+ }\r
+ }\r
+ else\r
+ {\r
+ oNorm = 0.;\r
+ oDot1 = 0.;\r
+ oDot2 = 0.;\r
+ }\r
+\r
+ oId = t_id;\r
+ return true;\r
+ }\r
+ }\r
+ ++nNum;\r
+ }\r
+\r
+ return false;\r
+}\r
+\r
+template <typename TInput, typename TOutput>\r
+bool\r
+FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::CheckTopology(OutputMeshType * oMesh, const NodeType & iNode)\r
+{\r
+ (void)oMesh;\r
+ (void)iNode;\r
+\r
+ return true;\r
+}\r
+\r
+template <typename TInput, typename TOutput>\r
+void\r
+FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::InitializeOutput(OutputMeshType * oMesh)\r
+{\r
+ printf("PG FastMarchingQuadEdgeMeshFilterResultsBase::InitializeOutput Start \n");\r
+ if(this->GetInput() != m_InputMesh || this->GetOutput() == NULL){\r
+ this->CopyInputMeshToOutputMeshGeometry();\r
+ }\r
+ else{\r
+ printf("PG FastMarchingQuadEdgeMeshFilterResultsBase::InitializeOutput FLAG 0.5 \n");\r
+ oMesh = this->GetOutput();\r
+ }\r
+ printf("PG FastMarchingQuadEdgeMeshFilterResultsBase::InitializeOutput FLAG 1 \n");\r
+ // Check that the input mesh is made of triangles\r
+ //We already now it is made of triangles\r
+ /**\r
+ {\r
+ OutputCellsContainerPointer cells = oMesh->GetCells();\r
+ OutputCellsContainerConstIterator c_it = cells->Begin();\r
+ OutputCellsContainerConstIterator c_end = cells->End();\r
+\r
+ while (c_it != c_end)\r
+ {\r
+ OutputCellType * cell = c_it.Value();\r
+\r
+ if (cell->GetNumberOfPoints() != 3)\r
+ {\r
+ itkGenericExceptionMacro(<< "Input mesh has non triangular faces");\r
+ }\r
+ ++c_it;\r
+ }\r
+ }*/\r
+\r
+ OutputPointsContainerPointer points = oMesh->GetPoints();\r
+\r
+ OutputPointDataContainerPointer pointdata = OutputPointDataContainer::New();\r
+ pointdata->Reserve(points->Size());\r
+\r
+ OutputPointsContainerIterator p_it = points->Begin();\r
+ OutputPointsContainerIterator p_end = points->End();\r
+\r
+ while (p_it != p_end)\r
+ {\r
+ pointdata->SetElement(p_it->Index(), this->m_LargeValue);\r
+ ++p_it;\r
+ }\r
+\r
+ oMesh->SetPointData(pointdata);\r
+\r
+ m_Label.clear();\r
+\r
+ if (this->m_AlivePoints)\r
+ {\r
+ NodePairContainerConstIterator pointsIter = this->m_AlivePoints->Begin();\r
+ NodePairContainerConstIterator pointsEnd = this->m_AlivePoints->End();\r
+ while (pointsIter != pointsEnd)\r
+ {\r
+ // get node from alive points container\r
+ NodeType idx = pointsIter->Value().GetNode();\r
+\r
+ if (points->IndexExists(idx))\r
+ {\r
+ OutputPixelType outputPixel = pointsIter->Value().GetValue();\r
+\r
+ this->SetLabelValueForGivenNode(idx, Traits::Alive);\r
+ this->SetOutputValue(oMesh, idx, outputPixel);\r
+ }\r
+\r
+ ++pointsIter;\r
+ }\r
+ }\r
+ if (this->m_ForbiddenPoints)\r
+ {\r
+ NodePairContainerConstIterator pointsIter = this->m_ForbiddenPoints->Begin();\r
+ NodePairContainerConstIterator pointsEnd = this->m_ForbiddenPoints->End();\r
+\r
+ OutputPixelType zero = NumericTraits<OutputPixelType>::ZeroValue();\r
+\r
+ while (pointsIter != pointsEnd)\r
+ {\r
+ NodeType idx = pointsIter->Value().GetNode();\r
+\r
+ if (points->IndexExists(idx))\r
+ {\r
+ this->SetLabelValueForGivenNode(idx, Traits::Forbidden);\r
+ this->SetOutputValue(oMesh, idx, zero);\r
+ }\r
+\r
+ ++pointsIter;\r
+ }\r
+ }\r
+ if (this->m_TrialPoints)\r
+ {\r
+ NodePairContainerConstIterator pointsIter = this->m_TrialPoints->Begin();\r
+ NodePairContainerConstIterator pointsEnd = this->m_TrialPoints->End();\r
+\r
+ while (pointsIter != pointsEnd)\r
+ {\r
+ NodeType idx = pointsIter->Value().GetNode();\r
+\r
+ if (points->IndexExists(idx))\r
+ {\r
+ OutputPixelType outputPixel = pointsIter->Value().GetValue();\r
+\r
+ this->SetLabelValueForGivenNode(idx, Traits::InitialTrial);\r
+ this->SetOutputValue(oMesh, idx, outputPixel);\r
+\r
+ this->m_Heap.push(pointsIter->Value());\r
+ }\r
+\r
+ ++pointsIter;\r
+ }\r
+ }printf("PG FastMarchingQuadEdgeMeshFilterResultsBase::InitializeOutput ENDED \n");\r
+}\r
+} // namespace itk\r
+\r
+#endif // itkFastMarchingQuadEdgeMeshFilterResultsBase_hxx\r