--- /dev/null
+//=====
+// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost)
+//=====
+#include "bbitkvtkGeodesicMeshDeformation.h"
+#include "bbitkvtkPackage.h"
+#include <vtkCellArrayIterator.h>
+namespace bbitkvtk
+{
+
+BBTK_ADD_BLACK_BOX_TO_PACKAGE(itkvtk,GeodesicMeshDeformation)
+BBTK_BLACK_BOX_IMPLEMENTATION(GeodesicMeshDeformation,bbtk::AtomicBlackBox);
+//=====
+// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost)
+//=====
+void GeodesicMeshDeformation::Process()
+{
+printf("PG GeodesicMeshDeformation::Process START \n");
+// THE MAIN PROCESSING METHOD BODY
+// Here we simply set the input 'In' value to the output 'Out'
+// And print out the output value
+// INPUT/OUTPUT ACCESSORS ARE OF THE FORM :
+// void bbSet{Input|Output}NAME(const TYPE&)
+// const TYPE& bbGet{Input|Output}NAME() const
+// Where :
+// * NAME is the name of the input/output
+// (the one provided in the attribute 'name' of the tag 'input')
+// * TYPE is the C++ type of the input/output
+// (the one provided in the attribute 'type' of the tag 'input')
+ std::vector<double> lstCenter = bbGetInputCenter();
+ double s = bbGetInputS();
+ bool ok = false;
+ using MeshType = itk::QuadEdgeMesh<double, 3>;
+
+ //Set up QuadEdge and filter every time polydata changes
+ if ((bbGetInputIn() != polydata) && (bbGetInputActive()==true) && (bbGetInputIn() != NULL))
+ {
+ //Reset displacement
+ if(lstCenter.size() != NULL){
+ backLstCenter[0] = lstCenter[0];
+ backLstCenter[1] = lstCenter[1];
+ backLstCenter[2] = lstCenter[2];
+ }
+ /**
+ Create QuadEdge
+ */
+ using PointType = MeshType::PointType;
+ //polydataBack
+ polydata = bbGetInputIn();
+
+ //points = points of input polydata
+ vtkPoints* points=bbGetInputIn()->GetPoints();
+
+ quadEdgeMesh = MeshType::New();
+ double currentPoint[3] = {};
+ PointType p;
+
+ int numberOfPoints = points->GetNumberOfPoints();;
+ for(int pointId = 0; pointId < numberOfPoints; pointId++){
+ points->GetPoint(pointId, currentPoint);
+ p[0] = currentPoint[0];
+ p[1] = currentPoint[1];
+ p[2] = currentPoint[2];
+ quadEdgeMesh->SetPoint(pointId, p);
+ //point data must be 1 on each point so that fast marching can calculate the distance.
+ quadEdgeMesh->SetPointData(pointId, 1.);
+ }
+ printf("PG GeodesicMeshDeformation::Process Entered FLAG 4 \n");
+ //add cells to QuadEdge mesh
+ auto cellIterator = vtk::TakeSmartPointer(bbGetInputIn()->GetPolys()->NewIterator());
+ vtkIdList* cell;
+ for(cellIterator->GoToFirstCell(); !cellIterator->IsDoneWithTraversal(); cellIterator->GoToNextCell()){
+ cell = cellIterator->GetCurrentCell();
+ quadEdgeMesh->AddFaceTriangle(cell->GetId(0), cell->GetId(1), cell->GetId(2));
+ }
+ printf("PG GeodesicMeshDeformation::Process FLAG 5\n");
+
+ using FastMarchingType = itk::FastMarchingQuadEdgeMeshFilterResultsBase<MeshType, MeshType>;
+ fmmFilter = FastMarchingType::New();
+ fmmFilter->SetInput(quadEdgeMesh);
+
+ using NodePairType = FastMarchingType::NodePairType;
+ using NodePairContainerType = FastMarchingType::NodePairContainerType;
+
+ auto trial = NodePairContainerType::New();
+ NodePairType nodePair(0, 0.);
+ trial->push_back(nodePair);
+
+ fmmFilter->SetTrialPoints(trial);
+
+ using CriterionType = itk::FastMarchingThresholdStoppingCriterion<MeshType, MeshType>;
+ auto criterion = CriterionType::New();
+ sCurrent = 1;
+ criterion->SetThreshold(1);
+ fmmFilter->SetStoppingCriterion(criterion);
+ fmmFilter->SetCollectPoints(true);
+ fmmFilter->SetReleaseDataBeforeUpdateFlag(false);
+ fmmFilter->Update();
+
+ printf("PG GeodesicMeshDeformation::Process Filter Created \n");
+ }
+
+ if (bbGetInputTypeIn()==0) // direction
+ {
+ if (bbGetInputDirection().size()==3)
+ {
+ ok = !( (bbGetInputDirection()[0]==0) && (bbGetInputDirection()[1]==0) && (bbGetInputDirection()[2]==0) );
+ }
+ } // if TypeIn 0
+
+ if (bbGetInputTypeIn()==1) // center
+ {
+ ok = ( lstCenter.size()==3 );
+ } // if TypeIn 1
+
+ if ( (bbGetInputIn()!=NULL) && (ok==true) && (bbGetInputEdgeId()>=0) && (bbGetInputActive()==true) )
+ {
+ vtkPoints* points=bbGetInputIn()->GetPoints();
+ long i,size=points->GetNumberOfPoints();
+ double p[3]; // point
+ double pb[3]; // point base
+ double np[3]; // new point
+ double sx,sy,sz;
+ sx = s*2;
+ sy = sx;
+ sz = sy;
+
+ double displcement_x = 0;
+ double displcement_y = 0;
+ double displcement_z = 0;
+
+ if (bbGetInputTypeIn()==0) // Direction
+ {
+ displcement_x = bbGetInputDirection()[0];
+ displcement_y = bbGetInputDirection()[1];
+ displcement_z = bbGetInputDirection()[2];
+ } // if TypeIn 0 Direction
+
+ printf(" EED GeodesicMeshDeformation::Process %ld %ld - %f %f %f \n", EdgeIdBack, bbGetInputEdgeId() , lstCenter[0],lstCenter[1],lstCenter[2] );
+
+ if (bbGetInputTypeIn()==1) // Center
+ {
+ if (EdgeIdBack==bbGetInputEdgeId() )
+ {
+ displcement_x = (lstCenter[0]-backLstCenter[0])/1.0;
+ displcement_y = (lstCenter[1]-backLstCenter[1])/1.0;
+ displcement_z = (lstCenter[2]-backLstCenter[2])/1.0;
+ } // if EdgeIdBack!=bbGetInputEdgeId()
+ backLstCenter[0] = lstCenter[0];
+ backLstCenter[1] = lstCenter[1];
+ backLstCenter[2] = lstCenter[2];
+
+ } // if TypeIn 1 Center
+ points->GetPoint( bbGetInputEdgeId() , pb );
+ if (EdgeIdBack!=bbGetInputEdgeId() )
+ {
+ EdgeIdBack = bbGetInputEdgeId();
+
+ backLstCenter[0] = lstCenter[0];
+ backLstCenter[1] = lstCenter[1];
+ backLstCenter[2] = lstCenter[2];
+
+ using NodePairType = FastMarchingType::NodePairType;
+ using NodePairContainerType = FastMarchingType::NodePairContainerType;
+
+ auto processedPoints = NodePairContainerType::New();
+ fmmFilter->SetProcessedPoints(processedPoints);
+
+ auto trial = NodePairContainerType::New();
+
+ //Seed point for fast marching with distance 0
+ NodePairType nodePair(bbGetInputEdgeId(), 0.);
+ trial->push_back(nodePair);
+
+ fmmFilter->SetTrialPoints(trial);
+
+ //threshold distance for fast marching method
+ if(sCurrent != s){
+ using CriterionType = itk::FastMarchingThresholdStoppingCriterion<MeshType, MeshType>;
+ auto criterion = CriterionType::New();
+ sCurrent = s;
+ criterion->SetThreshold(s*2);
+ fmmFilter->SetStoppingCriterion(criterion);
+ }
+ fmmFilter->Update();
+
+ } // if EdgeIdBack
+
+ if ( !((displcement_x==0) &&(displcement_y==0) && (displcement_z==0)) )
+ {
+ if(fmmFilter != NULL){
+ MeshType::PointDataContainer::Pointer pointData = fmmFilter->GetOutput()->GetPointData();
+ double pPD[3];
+
+ FastMarchingType::NodePairContainerType::Iterator BegProcessedIt = fmmFilter->GetProcessedPoints()->Begin();
+ FastMarchingType::NodePairContainerType::Iterator EndProcessedIt = fmmFilter->GetProcessedPoints()->End();
+
+ while (BegProcessedIt != EndProcessedIt)
+ {
+ //point modification
+ points->GetPoint(BegProcessedIt.Value().GetNode(), pPD);
+ double weight = exp(-BegProcessedIt.Value().GetValue() * 1. / s);
+ if(weight > 0.0001){
+ np[0] = pPD[0]+ weight*displcement_x;
+ np[1] = pPD[1]+ weight*displcement_y;
+ np[2] = pPD[2]+ weight*displcement_z;
+ points->SetPoint(BegProcessedIt.Value().GetNode(), np);
+ }
+ ++BegProcessedIt;
+ }
+ points->Modified();
+ bbGetInputIn()->Modified();
+ }
+ } // if distance_x y z != 0
+ } // In != NULL ok active
+ printf("PG GeodesicMeshDeformation::Process END \n");
+}
+//=====
+// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost)
+//=====
+void GeodesicMeshDeformation::bbUserSetDefaultValues()
+{
+
+// SET HERE THE DEFAULT INPUT/OUTPUT VALUES OF YOUR BOX
+// Here we initialize the input 'In' to 0
+ bbSetInputActive(true);
+ bbSetInputTypeIn(0);
+ bbSetInputIn(NULL);
+ std::vector<double> direction;
+ direction.push_back(1);
+ direction.push_back(0);
+ direction.push_back(0);
+ bbSetInputDirection(direction);
+ EdgeIdBack=-1;
+ bbSetInputEdgeId(EdgeIdBack);
+ bbSetInputS(10);
+
+ backLstCenter.push_back(0);
+ backLstCenter.push_back(0);
+ backLstCenter.push_back(0);
+ polydata = NULL;
+ quadEdgeMesh = NULL;
+ sCurrent = 0;
+ fmmFilter = NULL;
+ firstTime = true;
+}
+//=====
+// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost)
+//=====
+void GeodesicMeshDeformation::bbUserInitializeProcessing()
+{
+
+// THE INITIALIZATION METHOD BODY :
+// Here does nothing
+// but this is where you should allocate the internal/output pointers
+// if any
+
+
+}
+//=====
+// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost)
+//=====
+void GeodesicMeshDeformation::bbUserFinalizeProcessing()
+{
+
+// THE FINALIZATION METHOD BODY :
+// Here does nothing
+// but this is where you should desallocate the internal/output pointers
+// if any
+
+}
+}
+// EO namespace bbitkvtk
+
+
--- /dev/null
+//=====
+// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost)
+//=====
+#ifndef __bbitkvtkGeodesicMeshDeformation_h_INCLUDED__
+#define __bbitkvtkGeodesicMeshDeformation_h_INCLUDED__
+
+#include "bbitkvtk_EXPORT.h"
+#include "bbtkAtomicBlackBox.h"
+#include "iostream"
+
+#include "vtkPolyData.h"
+#include "itkQuadEdgeMesh.h"
+#include "itkFastMarchingThresholdStoppingCriterion.h"
+//#include "itkFastMarchingQuadEdgeMeshFilterBase.h"
+#include "itkFastMarchingQuadEdgeMeshFilterResultsBase.h"
+
+namespace bbitkvtk
+{
+
+class bbitkvtk_EXPORT GeodesicMeshDeformation
+ :
+ public bbtk::AtomicBlackBox
+{
+ BBTK_BLACK_BOX_INTERFACE(GeodesicMeshDeformation,bbtk::AtomicBlackBox);
+//=====
+// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost)
+//=====
+ BBTK_DECLARE_INPUT(Active, bool);
+ BBTK_DECLARE_INPUT(TypeIn, int);
+ BBTK_DECLARE_INPUT(In,vtkPolyData*);
+ BBTK_DECLARE_INPUT(EdgeId, long);
+ BBTK_DECLARE_INPUT(S, double);
+ BBTK_DECLARE_INPUT(Center, std::vector<double>);
+ BBTK_DECLARE_INPUT(Direction, std::vector<double>);
+ //BBTK_DECLARE_OUTPUT(Out,double);
+ BBTK_PROCESS(Process);
+ void Process();
+
+ using MeshType = itk::QuadEdgeMesh<double, 3>;
+ using FastMarchingType = itk::FastMarchingQuadEdgeMeshFilterResultsBase<MeshType, MeshType>;
+
+ long EdgeIdBack;
+ std::vector<int> voiIdPoints;
+ std::vector<double> backLstCenter;
+ vtkPolyData *polydata;
+ MeshType::Pointer quadEdgeMesh;
+ double sCurrent;
+ FastMarchingType::Pointer fmmFilter;
+ bool firstTime;
+//=====
+// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost)
+//=====
+};
+
+BBTK_BEGIN_DESCRIBE_BLACK_BOX(GeodesicMeshDeformation,bbtk::AtomicBlackBox);
+BBTK_NAME("GeodesicMeshDeformation");
+BBTK_AUTHOR("InfoDev");
+BBTK_DESCRIPTION("No Description.");
+BBTK_CATEGORY("empty");
+BBTK_INPUT(GeodesicMeshDeformation,Active,"(default true) true/false",bool,"");
+BBTK_INPUT(GeodesicMeshDeformation,TypeIn,"(default 0) 0:Direction 1:Center",int,"");
+BBTK_INPUT(GeodesicMeshDeformation,In,"vtk PolyData",vtkPolyData*,"");
+BBTK_INPUT(GeodesicMeshDeformation,EdgeId,"Edge Id",long,"");
+BBTK_INPUT(GeodesicMeshDeformation,S,"Deformation",double,"");
+BBTK_INPUT(GeodesicMeshDeformation,Center,"[X,Y,Z]",std::vector<double>,"");
+BBTK_INPUT(GeodesicMeshDeformation,Direction,"(default [1,0,0]) [X,Y,Z]",std::vector<double>,"");
+BBTK_END_DESCRIBE_BLACK_BOX(GeodesicMeshDeformation);
+//=====
+// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost)
+//=====
+}
+// EO namespace bbitkvtk
+
+#endif // __bbitkvtkGeodesicMeshDeformation_h_INCLUDED__
+
--- /dev/null
+/*=========================================================================\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 itkFastMarchingQuadEdgeMeshFilterBase_h\r
+#define itkFastMarchingQuadEdgeMeshFilterBase_h\r
+\r
+#include "itkFastMarchingBase.h"\r
+#include "itkFastMarchingTraits.h"\r
+\r
+namespace itk\r
+{\r
+/**\r
+ *\class FastMarchingQuadEdgeMeshFilterBase\r
+ \brief Fast Marching Method on QuadEdgeMesh\r
+\r
+ The speed function is specified by the input mesh. Data associated to each\r
+ point is considered as the speed function. The speed function is set using\r
+ the method SetInput().\r
+\r
+ If the speed function is constant and of value one, fast marching results is\r
+ an approximate geodesic function from the initial alive points.\r
+\r
+ Implementation of this class is based on\r
+ "Fast Marching Methods on Triangulated Domains", Kimmel, R., and Sethian, J.A.,\r
+ Proc. Nat. Acad. Sci., 95, pp. 8341-8435, 1998.\r
+\r
+ \ingroup ITKFastMarching\r
+*/\r
+template <typename TInput, typename TOutput>\r
+//2023-04-27-PG class ITK_TEMPLATE_EXPORT FastMarchingQuadEdgeMeshFilterBase : public FastMarchingBase<TInput, TOutput>\r
+class ITK_TEMPLATE_EXPORT FastMarchingQuadEdgeMeshFilterResultsBase : public FastMarchingBase<TInput, TOutput>\r
+{\r
+public:\r
+ //ITK_DISALLOW_COPY_AND_MOVE(FastMarchingQuadEdgeMeshFilterBase);\r
+ ITK_DISALLOW_COPY_AND_MOVE(FastMarchingQuadEdgeMeshFilterResultsBase);\r
+\r
+ //using Self = FastMarchingQuadEdgeMeshFilterBase;\r
+ using Self =FastMarchingQuadEdgeMeshFilterResultsBase;\r
+ using Superclass = FastMarchingBase<TInput, TOutput>;\r
+ using Pointer = SmartPointer<Self>;\r
+ using ConstPointer = SmartPointer<const Self>;\r
+ using Traits = typename Superclass::Traits;\r
+\r
+ /** Method for creation through the object factory. */\r
+ itkNewMacro(Self);\r
+\r
+ /** Run-time type information (and related methods). */\r
+ //itkTypeMacro(FastMarchingQuadEdgeMeshFilterBase, FastMarchingBase);\r
+ itkTypeMacro(FastMarchingQuadEdgeMeshFilterResultsBase, FastMarchingBase);\r
+\r
+ using InputMeshType = typename Superclass::InputDomainType;\r
+ using InputMeshPointer = typename Superclass::InputDomainPointer;\r
+ using InputPixelType = typename Superclass::InputPixelType;\r
+ using InputPointType = typename InputMeshType::PointType;\r
+ using InputPointIdentifierType = typename InputMeshType::PointIdentifier;\r
+\r
+ using OutputMeshType = typename Superclass::OutputDomainType;\r
+ using OutputMeshPointer = typename Superclass::OutputDomainPointer;\r
+ using OutputPixelType = typename Superclass::OutputPixelType;\r
+ using OutputPointType = typename OutputMeshType::PointType;\r
+ using OutputVectorType = typename OutputPointType::VectorType;\r
+ using OutputVectorRealType = typename OutputVectorType::RealValueType;\r
+ using OutputQEType = typename OutputMeshType::QEType;\r
+ using OutputPointIdentifierType = typename OutputMeshType::PointIdentifier;\r
+ using OutputPointsContainer = typename OutputMeshType::PointsContainer;\r
+ using OutputPointsContainerPointer = typename OutputPointsContainer::Pointer;\r
+ using OutputPointsContainerIterator = typename OutputPointsContainer::Iterator;\r
+ using OutputPointDataContainer = typename OutputMeshType::PointDataContainer;\r
+ using OutputPointDataContainerPointer = typename OutputPointDataContainer::Pointer;\r
+\r
+ using OutputCellsContainer = typename OutputMeshType::CellsContainer;\r
+ using OutputCellsContainerPointer = typename OutputCellsContainer::Pointer;\r
+ using OutputCellsContainerConstIterator = typename OutputCellsContainer::ConstIterator;\r
+ using OutputCellType = typename OutputMeshType::CellType;\r
+\r
+\r
+ using NodeType = typename Traits::NodeType;\r
+ using NodePairType = typename Traits::NodePairType;\r
+ using NodePairContainerType = typename Traits::NodePairContainerType;\r
+ using NodePairContainerPointer = typename Traits::NodePairContainerPointer;\r
+ using NodePairContainerConstIterator = typename Traits::NodePairContainerConstIterator;\r
+\r
+ // using NodeContainerType = typename Traits::NodeContainerType;\r
+ // using NodeContainerPointer = typename Traits::NodeContainerPointer;\r
+ // using NodeContainerConstIterator = typename Traits::NodeContainerConstIterator;\r
+\r
+ using LabelType = typename Superclass::LabelType;\r
+\r
+ using NodeLabelMapType = std::map<NodeType, LabelType>;\r
+ using NodeLabelMapIterator = typename NodeLabelMapType::iterator;\r
+ using NodeLabelMapConstIterator = typename NodeLabelMapType::const_iterator;\r
+\r
+protected:\r
+ //FastMarchingQuadEdgeMeshFilterBase();\r
+ FastMarchingQuadEdgeMeshFilterResultsBase();\r
+ //~FastMarchingQuadEdgeMeshFilterBase() override = default;\r
+ ~FastMarchingQuadEdgeMeshFilterResultsBase() override = default;\r
+\r
+ NodeLabelMapType m_Label;\r
+\r
+ IdentifierType\r
+ GetTotalNumberOfNodes() const override;\r
+\r
+ void\r
+ SetOutputValue(OutputMeshType * oMesh, const NodeType & iNode, const OutputPixelType & iValue) override;\r
+\r
+ const OutputPixelType\r
+ GetOutputValue(OutputMeshType * oMesh, const NodeType & iNode) const override;\r
+\r
+ unsigned char\r
+ GetLabelValueForGivenNode(const NodeType & iNode) const override;\r
+\r
+ void\r
+ SetLabelValueForGivenNode(const NodeType & iNode, const LabelType & iLabel) override;\r
+\r
+ void\r
+ UpdateNeighbors(OutputMeshType * oMesh, const NodeType & iNode) override;\r
+\r
+ void\r
+ UpdateValue(OutputMeshType * oMesh, const NodeType & iNode) override;\r
+\r
+ const OutputVectorRealType\r
+ 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
+\r
+ const OutputVectorRealType\r
+ 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
+ bool\r
+ 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
+ bool\r
+ CheckTopology(OutputMeshType * oMesh, const NodeType & iNode) override;\r
+\r
+ void\r
+ InitializeOutput(OutputMeshType * oMesh) override;\r
+\r
+private:\r
+ const InputMeshType * m_InputMesh;\r
+ const InputMeshType * m_BackInputMesh;\r
+};\r
+} // namespace itk\r
+\r
+#ifndef ITK_MANUAL_INSTANTIATION\r
+//# include "itkFastMarchingQuadEdgeMeshFilterBase.hxx"\r
+# include "itkFastMarchingQuadEdgeMeshFilterResultsBase.hxx"\r
+#endif\r
+\r
+#endif // itkFastMarchingQuadEdgeMeshFilterBase_h\r
--- /dev/null
+/*=========================================================================\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
git://git.creatis.insa-lyon.fr / bbtk.git / commitdiff