packages/itkvtk/src/itkFastMarchingQuadEdgeMeshFilterResultsBase.hxx

   1 /*=========================================================================\r
   2  *\r
   3  *  Copyright NumFOCUS\r
   4  *\r
   5  *  Licensed under the Apache License, Version 2.0 (the "License");\r
   6  *  you may not use this file except in compliance with the License.\r
   7  *  You may obtain a copy of the License at\r
   8  *\r
   9  *         http://www.apache.org/licenses/LICENSE-2.0.txt\r
  10  *\r
  11  *  Unless required by applicable law or agreed to in writing, software\r
  12  *  distributed under the License is distributed on an "AS IS" BASIS,\r
  13  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\r
  14  *  See the License for the specific language governing permissions and\r
  15  *  limitations under the License.\r
  16  *\r
  17  *=========================================================================*/\r
  18 \r
  19 #ifndef itkFastMarchingQuadEdgeMeshFilterResultsBase_hxx\r
  20 #define itkFastMarchingQuadEdgeMeshFilterResultsBase_hxx\r
  21 \r
  22 #include "itkFastMarchingQuadEdgeMeshFilterResultsBase.h"\r
  23 \r
  24 #include "itkMath.h"\r
  25 #include "itkVector.h"\r
  26 #include "itkMatrix.h"\r
  27 \r
  28 namespace itk\r
  29 {\r
  30 template <typename TInput, typename TOutput>\r
  31 FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::FastMarchingQuadEdgeMeshFilterResultsBase()\r
  32   : Superclass()\r
  33 {\r
  34   this->m_InputMesh = nullptr;\r
  35   this->m_BackInputMesh = nullptr;\r
  36 }\r
  37 \r
  38 template <typename TInput, typename TOutput>\r
  39 IdentifierType\r
  40 FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::GetTotalNumberOfNodes() const\r
  41 {\r
  42   return this->GetInput()->GetNumberOfPoints();\r
  43 }\r
  44 \r
  45 template <typename TInput, typename TOutput>\r
  46 void\r
  47 FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::SetOutputValue(OutputMeshType *        oMesh,\r
  48                                                                     const NodeType &        iNode,\r
  49                                                                     const OutputPixelType & iValue)\r
  50 {\r
  51   oMesh->SetPointData(iNode, iValue);\r
  52 }\r
  53 \r
  54 template <typename TInput, typename TOutput>\r
  55 const typename FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::OutputPixelType\r
  56 FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::GetOutputValue(OutputMeshType * oMesh,\r
  57                                                                     const NodeType & iNode) const\r
  58 {\r
  59   OutputPixelType outputValue = NumericTraits<OutputPixelType>::ZeroValue();\r
  60   oMesh->GetPointData(iNode, &outputValue);\r
  61   return outputValue;\r
  62 }\r
  63 \r
  64 \r
  65 template <typename TInput, typename TOutput>\r
  66 unsigned char\r
  67 FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::GetLabelValueForGivenNode(const NodeType & iNode) const\r
  68 {\r
  69   auto it = m_Label.find(iNode);\r
  70 \r
  71   if (it != m_Label.end())\r
  72   {\r
  73     return it->second;\r
  74   }\r
  75   else\r
  76   {\r
  77     return Traits::Far;\r
  78   }\r
  79 }\r
  80 \r
  81 template <typename TInput, typename TOutput>\r
  82 void\r
  83 FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::SetLabelValueForGivenNode(const NodeType &  iNode,\r
  84                                                                                const LabelType & iLabel)\r
  85 {\r
  86   m_Label[iNode] = iLabel;\r
  87 }\r
  88 \r
  89 template <typename TInput, typename TOutput>\r
  90 void\r
  91 FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::UpdateNeighbors(OutputMeshType * oMesh, const NodeType & iNode)\r
  92 {\r
  93   //printf("PG FastMarchingQuadEdgeMeshFilterResultsBase::UpdateNeighbors Start \n");\r
  94   OutputPointType p;\r
  95   oMesh->GetPoint(iNode, &p);\r
  96 \r
  97   OutputQEType * qe = p.GetEdge();\r
  98 \r
  99   if (qe)\r
 100   {\r
 101     OutputQEType * qe_it = qe;\r
 102     this->m_InputMesh = this->GetInput();\r
 103     do\r
 104     {\r
 105       if (qe_it)\r
 106       {\r
 107         OutputPointIdentifierType neigh_id = qe_it->GetDestination();\r
 108 \r
 109         const char label = this->GetLabelValueForGivenNode(neigh_id);\r
 110 \r
 111         if ((label != Traits::Alive) && (label != Traits::InitialTrial) && (label != Traits::Forbidden))\r
 112         {\r
 113           this->UpdateValue(oMesh, neigh_id);\r
 114         }\r
 115       }\r
 116       else\r
 117       {\r
 118         itkGenericExceptionMacro(<< "qe_it is nullptr");\r
 119       }\r
 120       qe_it = qe_it->GetOnext();\r
 121     } while (qe_it != qe);\r
 122   }\r
 123   else\r
 124   {\r
 125     itkGenericExceptionMacro(<< "qe is nullptr");\r
 126   }\r
 127   //printf("PG FastMarchingQuadEdgeMeshFilterResultsBase::UpdateNeighbors End \n");\r
 128 }\r
 129 \r
 130 template <typename TInput, typename TOutput>\r
 131 void\r
 132 FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::UpdateValue(OutputMeshType * oMesh, const NodeType & iNode)\r
 133 {\r
 134   OutputPointType p;\r
 135   oMesh->GetPoint(iNode, &p);\r
 136 \r
 137   InputPixelType F = NumericTraits<InputPixelType>::ZeroValue();\r
 138   //this->m_InputMesh->GetPointData(iNode, &F);\r
 139   F = 1.;\r
 140   if (F < 0.)\r
 141   {\r
 142     itkGenericExceptionMacro(<< "F < 0.");\r
 143   }\r
 144 \r
 145   OutputQEType * qe = p.GetEdge();\r
 146 \r
 147   OutputPixelType outputPixel = this->m_LargeValue;\r
 148 \r
 149   if (qe)\r
 150   {\r
 151     OutputQEType * qe_it = qe;\r
 152     qe_it = qe_it->GetOnext();\r
 153 \r
 154     do\r
 155     {\r
 156       OutputQEType * qe_it2 = qe_it->GetOnext();\r
 157 \r
 158       if (qe_it2)\r
 159       {\r
 160         if (qe_it->GetLeft() != OutputMeshType::m_NoFace)\r
 161         {\r
 162           OutputPointIdentifierType id1 = qe_it->GetDestination();\r
 163           OutputPointIdentifierType id2 = qe_it2->GetDestination();\r
 164 \r
 165           const auto label1 = static_cast<LabelType>(this->GetLabelValueForGivenNode(id1));\r
 166           const auto label2 = static_cast<LabelType>(this->GetLabelValueForGivenNode(id2));\r
 167 \r
 168           bool IsFar1 = (label1 != Traits::Far);\r
 169           bool IsFar2 = (label2 != Traits::Far);\r
 170 \r
 171           if (IsFar1 || IsFar2)\r
 172           {\r
 173             OutputPointType q1 = oMesh->GetPoint(id1);\r
 174             OutputPointType q2 = oMesh->GetPoint(id2);\r
 175 \r
 176             auto val1 = static_cast<OutputVectorRealType>(this->GetOutputValue(oMesh, id1));\r
 177             auto val2 = static_cast<OutputVectorRealType>(this->GetOutputValue(oMesh, id2));\r
 178 \r
 179             if (val1 > val2)\r
 180             {\r
 181               OutputPointType temp_pt(q1);\r
 182               q1 = q2;\r
 183               q2 = temp_pt;\r
 184 \r
 185               std::swap(val1, val2);\r
 186               std::swap(IsFar1, IsFar2);\r
 187               std::swap(id1, id2);\r
 188             }\r
 189 \r
 190             const OutputVectorRealType temp =\r
 191               this->Solve(oMesh, iNode, p, F, id1, q1, IsFar1, val1, id2, q2, IsFar2, val2);\r
 192 \r
 193             outputPixel = std::min(outputPixel, static_cast<OutputPixelType>(temp));\r
 194           }\r
 195         }\r
 196 \r
 197         qe_it = qe_it2;\r
 198       }\r
 199       else\r
 200       {\r
 201         // throw one exception here\r
 202         itkGenericExceptionMacro(<< "qe_it2 is nullptr");\r
 203       }\r
 204     } while (qe_it != qe);\r
 205 \r
 206     if (outputPixel < this->m_LargeValue)\r
 207     {\r
 208       this->SetOutputValue(oMesh, iNode, outputPixel);\r
 209 \r
 210       this->SetLabelValueForGivenNode(iNode, Traits::Trial);\r
 211 \r
 212       this->m_Heap.push(NodePairType(iNode, outputPixel));\r
 213     }\r
 214   }\r
 215   else\r
 216   {\r
 217     // throw one exception\r
 218     itkGenericExceptionMacro(<< "qe_it is nullptr");\r
 219   }\r
 220 }\r
 221 \r
 222 template <typename TInput, typename TOutput>\r
 223 const typename FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::OutputVectorRealType\r
 224 FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::Solve(OutputMeshType *             oMesh,\r
 225                                                            const NodeType &             iId,\r
 226                                                            const OutputPointType &      iCurrentPoint,\r
 227                                                            const OutputVectorRealType & iF,\r
 228                                                            const NodeType &             iId1,\r
 229                                                            const OutputPointType &      iP1,\r
 230                                                            const bool &                 iIsFar1,\r
 231                                                            const OutputVectorRealType   iVal1,\r
 232                                                            const NodeType &             iId2,\r
 233                                                            const OutputPointType &      iP2,\r
 234                                                            const bool &                 iIsFar2,\r
 235                                                            const OutputVectorRealType & iVal2) const\r
 236 {\r
 237   OutputVectorType Edge1 = iP1 - iCurrentPoint;\r
 238   OutputVectorType Edge2 = iP2 - iCurrentPoint;\r
 239 \r
 240   OutputVectorRealType sq_norm1 = Edge1.GetSquaredNorm();\r
 241   OutputVectorRealType norm1 = 0.;\r
 242 \r
 243   OutputVectorRealType epsilon = NumericTraits<OutputVectorRealType>::epsilon();\r
 244 \r
 245   if (sq_norm1 > epsilon)\r
 246   {\r
 247     norm1 = std::sqrt(sq_norm1);\r
 248 \r
 249     OutputVectorRealType inv_norm1 = 1. / norm1;\r
 250     Edge1 *= inv_norm1;\r
 251   }\r
 252 \r
 253   OutputVectorRealType sq_norm2 = Edge2.GetSquaredNorm();\r
 254   OutputVectorRealType norm2 = 0.;\r
 255 \r
 256   if (sq_norm2 > epsilon)\r
 257   {\r
 258     norm2 = std::sqrt(sq_norm2);\r
 259 \r
 260     OutputVectorRealType inv_norm2 = 1. / norm2;\r
 261     Edge2 *= inv_norm2;\r
 262   }\r
 263 \r
 264   if (!iIsFar1 && iIsFar2)\r
 265   {\r
 266     // only one point is a contributor\r
 267     return iVal2 + norm2 * iF;\r
 268   }\r
 269   if (iIsFar1 && !iIsFar2)\r
 270   {\r
 271     // only one point is a contributor\r
 272     return iVal1 + norm1 * iF;\r
 273   }\r
 274 \r
 275   if (iIsFar1 && iIsFar2)\r
 276   {\r
 277     auto dot = static_cast<OutputVectorRealType>(Edge1 * Edge2);\r
 278 \r
 279     if (dot >= 0.)\r
 280     {\r
 281       return ComputeUpdate(iVal1, iVal2, norm2, sq_norm2, norm1, sq_norm1, dot, iF);\r
 282     }\r
 283     else\r
 284     {\r
 285       OutputVectorRealType      sq_norm3, norm3, dot1, dot2;\r
 286       OutputPointIdentifierType new_id;\r
 287 \r
 288       bool unfolded =\r
 289         UnfoldTriangle(oMesh, iId, iCurrentPoint, iId1, iP1, iId2, iP2, norm3, sq_norm3, dot1, dot2, new_id);\r
 290 \r
 291       if (unfolded)\r
 292       {\r
 293         OutputVectorRealType t_sq_norm3 = norm3 * norm3;\r
 294 \r
 295         auto                 val3 = static_cast<OutputVectorRealType>(this->GetOutputValue(oMesh, new_id));\r
 296         OutputVectorRealType t1 = ComputeUpdate(iVal1, val3, norm3, t_sq_norm3, norm1, sq_norm1, dot1, iF);\r
 297         OutputVectorRealType t2 = ComputeUpdate(iVal2, val3, norm3, t_sq_norm3, norm2, sq_norm2, dot2, iF);\r
 298 \r
 299         return std::min(t1, t2);\r
 300       }\r
 301       else\r
 302       {\r
 303         return ComputeUpdate(iVal1, iVal2, norm2, sq_norm2, norm1, sq_norm1, dot, iF);\r
 304       }\r
 305     }\r
 306   }\r
 307 \r
 308   return this->m_LargeValue;\r
 309 }\r
 310 \r
 311 template <typename TInput, typename TOutput>\r
 312 const typename FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::OutputVectorRealType\r
 313 FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::ComputeUpdate(const OutputVectorRealType & iVal1,\r
 314                                                                    const OutputVectorRealType & iVal2,\r
 315                                                                    const OutputVectorRealType & iNorm1,\r
 316                                                                    const OutputVectorRealType & iSqNorm1,\r
 317                                                                    const OutputVectorRealType & iNorm2,\r
 318                                                                    const OutputVectorRealType & iSqNorm2,\r
 319                                                                    const OutputVectorRealType & iDot,\r
 320                                                                    const OutputVectorRealType & iF) const\r
 321 {\r
 322   auto large_value = static_cast<OutputVectorRealType>(this->m_LargeValue);\r
 323 \r
 324   OutputVectorRealType t = large_value;\r
 325 \r
 326   OutputVectorRealType CosAngle = iDot;\r
 327   OutputVectorRealType SinAngle = std::sqrt(1. - iDot * iDot);\r
 328 \r
 329   OutputVectorRealType u = iVal2 - iVal1;\r
 330 \r
 331   OutputVectorRealType sq_u = u * u;\r
 332   OutputVectorRealType f2 = iSqNorm1 + iSqNorm2 - 2. * iNorm1 * iNorm2 * CosAngle;\r
 333   OutputVectorRealType f1 = iNorm2 * u * (iNorm1 * CosAngle - iNorm2);\r
 334   OutputVectorRealType f0 = iSqNorm2 * (sq_u - iF * iF * iSqNorm1 * SinAngle * SinAngle);\r
 335 \r
 336   OutputVectorRealType delta = f1 * f1 - f0 * f2;\r
 337 \r
 338   OutputVectorRealType epsilon = NumericTraits<OutputVectorRealType>::epsilon();\r
 339 \r
 340   if (delta >= 0.)\r
 341   {\r
 342     if (itk::Math::abs(f2) > epsilon)\r
 343     {\r
 344       t = (-f1 - std::sqrt(delta)) / f2;\r
 345 \r
 346       // test if we must must choose the other solution\r
 347       if ((t < u) || (iNorm2 * (t - u) / t < iNorm1 * CosAngle) || (iNorm1 / CosAngle < iNorm2 * (t - u) / t))\r
 348       {\r
 349         t = (-f1 + std::sqrt(delta)) / f2;\r
 350       }\r
 351     }\r
 352     else\r
 353     {\r
 354       if (itk::Math::abs(f1) > epsilon)\r
 355       {\r
 356         t = -f0 / f1;\r
 357       }\r
 358       else\r
 359       {\r
 360         t = -large_value;\r
 361       }\r
 362     }\r
 363   }\r
 364   else\r
 365   {\r
 366     t = -large_value;\r
 367   }\r
 368 \r
 369   // choose the update from the 2 vertex only if upwind criterion is met\r
 370   if ((u < t) && (iNorm1 * CosAngle < iNorm2 * (t - u) / t) && (iNorm2 * (t - u) / t < iNorm1 / CosAngle))\r
 371   {\r
 372     return t + iVal1;\r
 373   }\r
 374   else\r
 375   {\r
 376     return std::min(iNorm2 * iF + iVal1, iNorm1 * iF + iVal2);\r
 377   }\r
 378 }\r
 379 \r
 380 template <typename TInput, typename TOutput>\r
 381 bool\r
 382 FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::UnfoldTriangle(OutputMeshType *                  oMesh,\r
 383                                                                     const OutputPointIdentifierType & iId,\r
 384                                                                     const OutputPointType &           iP,\r
 385                                                                     const OutputPointIdentifierType & iId1,\r
 386                                                                     const OutputPointType &           iP1,\r
 387                                                                     const OutputPointIdentifierType & iId2,\r
 388                                                                     const OutputPointType &           iP2,\r
 389                                                                     OutputVectorRealType &            oNorm,\r
 390                                                                     OutputVectorRealType &            oSqNorm,\r
 391                                                                     OutputVectorRealType &            oDot1,\r
 392                                                                     OutputVectorRealType &            oDot2,\r
 393                                                                     OutputPointIdentifierType &       oId) const\r
 394 {\r
 395   (void)iId;\r
 396 \r
 397   OutputVectorType     Edge1 = iP1 - iP;\r
 398   OutputVectorRealType Norm1 = Edge1.GetNorm();\r
 399 \r
 400   OutputVectorRealType epsilon = NumericTraits<OutputVectorRealType>::epsilon();\r
 401 \r
 402   if (Norm1 > epsilon)\r
 403   {\r
 404     OutputVectorRealType inv_Norm = 1. / Norm1;\r
 405     Edge1 *= inv_Norm;\r
 406   }\r
 407 \r
 408   OutputVectorType     Edge2 = iP2 - iP;\r
 409   OutputVectorRealType Norm2 = Edge2.GetNorm();\r
 410   if (Norm2 > epsilon)\r
 411   {\r
 412     OutputVectorRealType inv_Norm = 1. / Norm2;\r
 413     Edge2 *= inv_Norm;\r
 414   }\r
 415 \r
 416   auto dot = static_cast<OutputVectorRealType>(Edge1 * Edge2);\r
 417 \r
 418   // the equation of the lines defining the unfolding region\r
 419   // [e.g. line 1 : {x ; <x,eq1>=0} ]\r
 420   using Vector2DType = Vector<OutputVectorRealType, 2>;\r
 421   using Matrix2DType = Matrix<OutputVectorRealType, 2, 2>;\r
 422 \r
 423   Vector2DType v1;\r
 424   v1[0] = dot;\r
 425   v1[1] = std::sqrt(1. - dot * dot);\r
 426 \r
 427   Vector2DType v2;\r
 428   v2[0] = 1.;\r
 429   v2[1] = 0.;\r
 430 \r
 431   Vector2DType x1;\r
 432   x1[0] = Norm1;\r
 433   x1[1] = 0.;\r
 434 \r
 435   Vector2DType x2 = Norm2 * v1;\r
 436 \r
 437   // keep track of the starting point\r
 438   Vector2DType x_start1(x1);\r
 439   Vector2DType x_start2(x2);\r
 440 \r
 441   OutputPointIdentifierType id1 = iId1;\r
 442   OutputPointIdentifierType id2 = iId2;\r
 443 \r
 444   OutputQEType * qe = oMesh->FindEdge(id1, id2);\r
 445   qe = qe->GetSym();\r
 446 \r
 447   OutputPointType t_pt1 = iP1;\r
 448   OutputPointType t_pt2 = iP2;\r
 449   OutputPointType t_pt = t_pt1;\r
 450 \r
 451   unsigned int nNum = 0;\r
 452   while (nNum < 50 && qe->GetLeft() != OutputMeshType::m_NoFace)\r
 453   {\r
 454     OutputQEType *            qe_Lnext = qe->GetLnext();\r
 455     OutputPointIdentifierType t_id = qe_Lnext->GetDestination();\r
 456 \r
 457     oMesh->GetPoint(t_id, &t_pt);\r
 458 \r
 459     Edge1 = t_pt2 - t_pt1;\r
 460     Norm1 = Edge1.GetNorm();\r
 461     if (Norm1 > epsilon)\r
 462     {\r
 463       OutputVectorRealType inv_Norm = 1. / Norm1;\r
 464       Edge1 *= inv_Norm;\r
 465     }\r
 466 \r
 467     Edge2 = t_pt - t_pt1;\r
 468     Norm2 = Edge2.GetNorm();\r
 469     if (Norm2 > epsilon)\r
 470     {\r
 471       OutputVectorRealType inv_Norm = 1. / Norm2;\r
 472       Edge2 *= inv_Norm;\r
 473     }\r
 474 \r
 475     /* compute the position of the new point x on the unfolding plane (via a rotation of -alpha on (x2-x1)/rNorm1 )\r
 476             | cos(alpha) sin(alpha)|\r
 477         x = |-sin(alpha) cos(alpha)| * [x2-x1]*rNorm2/rNorm1 + x1   where cos(alpha)=dot\r
 478     */\r
 479     Vector2DType vv;\r
 480     if (Norm1 > epsilon)\r
 481     {\r
 482       vv = (x2 - x1) * Norm2 / Norm1;\r
 483     }\r
 484     else\r
 485     {\r
 486       vv.Fill(0.);\r
 487     }\r
 488 \r
 489     dot = Edge1 * Edge2;\r
 490 \r
 491     Matrix2DType rotation;\r
 492     rotation[0][0] = dot;\r
 493     rotation[0][1] = std::sqrt(1. - dot * dot);\r
 494     rotation[1][0] = -rotation[0][1];\r
 495     rotation[1][1] = dot;\r
 496 \r
 497     Vector2DType x = rotation * vv + x1;\r
 498 \r
 499 \r
 500     /* compute the intersection points.\r
 501        We look for x=x1+lambda*(x-x1) or x=x2+lambda*(x-x2) with <x,eqi>=0, so */\r
 502     OutputVectorRealType lambda11 = -(x1 * v1) / ((x - x1) * v1); // left most\r
 503     OutputVectorRealType lambda12 = -(x1 * v2) / ((x - x1) * v2); // right most\r
 504     OutputVectorRealType lambda21 = -(x2 * v1) / ((x - x2) * v1); // left most\r
 505     OutputVectorRealType lambda22 = -(x2 * v2) / ((x - x2) * v2); // right most\r
 506     bool                 bIntersect11 = (lambda11 >= 0.) && (lambda11 <= 1.);\r
 507     bool                 bIntersect12 = (lambda12 >= 0.) && (lambda12 <= 1.);\r
 508     bool                 bIntersect21 = (lambda21 >= 0.) && (lambda21 <= 1.);\r
 509     bool                 bIntersect22 = (lambda22 >= 0.) && (lambda22 <= 1.);\r
 510 \r
 511     if (bIntersect11 && bIntersect12)\r
 512     {\r
 513       qe = oMesh->FindEdge(id1, t_id);\r
 514       qe = qe->GetSym();\r
 515       id2 = t_id;\r
 516       t_pt2 = t_pt;\r
 517       x2 = x;\r
 518     }\r
 519     else\r
 520     {\r
 521       if (bIntersect21 && bIntersect22)\r
 522       {\r
 523         qe = oMesh->FindEdge(id2, t_id);\r
 524         qe = qe->GetSym();\r
 525         id1 = t_id;\r
 526         t_pt1 = t_pt;\r
 527         x1 = x;\r
 528       }\r
 529       else\r
 530       { // bIntersect11 && !bIntersect12 && !bIntersect21 && bIntersect22\r
 531 \r
 532         // that's it, we have found the point\r
 533         oSqNorm = x.GetSquaredNorm();\r
 534 \r
 535         if (oSqNorm > epsilon)\r
 536         {\r
 537           oNorm = std::sqrt(oSqNorm);\r
 538           OutputVectorRealType temp_norm = x_start1.GetNorm();\r
 539           if (temp_norm > epsilon)\r
 540           {\r
 541             oDot1 = x * x_start1 / (oNorm * temp_norm);\r
 542           }\r
 543           else\r
 544           {\r
 545             oDot1 = 0.;\r
 546           }\r
 547           temp_norm = x_start2.GetNorm();\r
 548           if (temp_norm > epsilon)\r
 549           {\r
 550             oDot2 = x * x_start2 / (oNorm * temp_norm);\r
 551           }\r
 552           else\r
 553           {\r
 554             oDot2 = 0.;\r
 555           }\r
 556         }\r
 557         else\r
 558         {\r
 559           oNorm = 0.;\r
 560           oDot1 = 0.;\r
 561           oDot2 = 0.;\r
 562         }\r
 563 \r
 564         oId = t_id;\r
 565         return true;\r
 566       }\r
 567     }\r
 568     ++nNum;\r
 569   }\r
 570 \r
 571   return false;\r
 572 }\r
 573 \r
 574 template <typename TInput, typename TOutput>\r
 575 bool\r
 576 FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::CheckTopology(OutputMeshType * oMesh, const NodeType & iNode)\r
 577 {\r
 578   (void)oMesh;\r
 579   (void)iNode;\r
 580 \r
 581   return true;\r
 582 }\r
 583 \r
 584 template <typename TInput, typename TOutput>\r
 585 void\r
 586 FastMarchingQuadEdgeMeshFilterResultsBase<TInput, TOutput>::InitializeOutput(OutputMeshType * oMesh)\r
 587 {\r
 588   printf("PG FastMarchingQuadEdgeMeshFilterResultsBase::InitializeOutput Start \n");\r
 589   if(this->GetInput() != m_InputMesh || this->GetOutput() == NULL){\r
 590           this->CopyInputMeshToOutputMeshGeometry();\r
 591   }\r
 592   else{\r
 593                 printf("PG FastMarchingQuadEdgeMeshFilterResultsBase::InitializeOutput FLAG 0.5 \n");\r
 594           oMesh = this->GetOutput();\r
 595   }\r
 596   printf("PG FastMarchingQuadEdgeMeshFilterResultsBase::InitializeOutput FLAG 1 \n");\r
 597   // Check that the input mesh is made of triangles\r
 598   //We already now it is made of triangles\r
 599   /**\r
 600   {\r
 601     OutputCellsContainerPointer       cells = oMesh->GetCells();\r
 602     OutputCellsContainerConstIterator c_it = cells->Begin();\r
 603     OutputCellsContainerConstIterator c_end = cells->End();\r
 604 \r
 605     while (c_it != c_end)\r
 606     {\r
 607       OutputCellType * cell = c_it.Value();\r
 608 \r
 609       if (cell->GetNumberOfPoints() != 3)\r
 610       {\r
 611         itkGenericExceptionMacro(<< "Input mesh has non triangular faces");\r
 612       }\r
 613       ++c_it;\r
 614     }\r
 615   }*/\r
 616 \r
 617   OutputPointsContainerPointer points = oMesh->GetPoints();\r
 618 \r
 619   OutputPointDataContainerPointer pointdata = OutputPointDataContainer::New();\r
 620   pointdata->Reserve(points->Size());\r
 621 \r
 622   OutputPointsContainerIterator p_it = points->Begin();\r
 623   OutputPointsContainerIterator p_end = points->End();\r
 624 \r
 625   while (p_it != p_end)\r
 626   {\r
 627     pointdata->SetElement(p_it->Index(), this->m_LargeValue);\r
 628     ++p_it;\r
 629   }\r
 630 \r
 631   oMesh->SetPointData(pointdata);\r
 632 \r
 633   m_Label.clear();\r
 634 \r
 635   if (this->m_AlivePoints)\r
 636   {\r
 637     NodePairContainerConstIterator pointsIter = this->m_AlivePoints->Begin();\r
 638     NodePairContainerConstIterator pointsEnd = this->m_AlivePoints->End();\r
 639     while (pointsIter != pointsEnd)\r
 640     {\r
 641       // get node from alive points container\r
 642       NodeType idx = pointsIter->Value().GetNode();\r
 643 \r
 644       if (points->IndexExists(idx))\r
 645       {\r
 646         OutputPixelType outputPixel = pointsIter->Value().GetValue();\r
 647 \r
 648         this->SetLabelValueForGivenNode(idx, Traits::Alive);\r
 649         this->SetOutputValue(oMesh, idx, outputPixel);\r
 650       }\r
 651 \r
 652       ++pointsIter;\r
 653     }\r
 654   }\r
 655   if (this->m_ForbiddenPoints)\r
 656   {\r
 657     NodePairContainerConstIterator pointsIter = this->m_ForbiddenPoints->Begin();\r
 658     NodePairContainerConstIterator pointsEnd = this->m_ForbiddenPoints->End();\r
 659 \r
 660     OutputPixelType zero = NumericTraits<OutputPixelType>::ZeroValue();\r
 661 \r
 662     while (pointsIter != pointsEnd)\r
 663     {\r
 664       NodeType idx = pointsIter->Value().GetNode();\r
 665 \r
 666       if (points->IndexExists(idx))\r
 667       {\r
 668         this->SetLabelValueForGivenNode(idx, Traits::Forbidden);\r
 669         this->SetOutputValue(oMesh, idx, zero);\r
 670       }\r
 671 \r
 672       ++pointsIter;\r
 673     }\r
 674   }\r
 675   if (this->m_TrialPoints)\r
 676   {\r
 677     NodePairContainerConstIterator pointsIter = this->m_TrialPoints->Begin();\r
 678     NodePairContainerConstIterator pointsEnd = this->m_TrialPoints->End();\r
 679 \r
 680     while (pointsIter != pointsEnd)\r
 681     {\r
 682       NodeType idx = pointsIter->Value().GetNode();\r
 683 \r
 684       if (points->IndexExists(idx))\r
 685       {\r
 686         OutputPixelType outputPixel = pointsIter->Value().GetValue();\r
 687 \r
 688         this->SetLabelValueForGivenNode(idx, Traits::InitialTrial);\r
 689         this->SetOutputValue(oMesh, idx, outputPixel);\r
 690 \r
 691         this->m_Heap.push(pointsIter->Value());\r
 692       }\r
 693 \r
 694       ++pointsIter;\r
 695     }\r
 696   }printf("PG FastMarchingQuadEdgeMeshFilterResultsBase::InitializeOutput ENDED \n");\r
 697 }\r
 698 } // namespace itk\r
 699 \r
 700 #endif // itkFastMarchingQuadEdgeMeshFilterResultsBase_hxx\r