From 4950aae691e2b5db557c397e5a34c49524c58e5a Mon Sep 17 00:00:00 2001 From: Pablo Garzon Date: Tue, 2 May 2023 16:17:11 +0200 Subject: [PATCH] #3501 Geodesic Deformation --- .../src/bbitkvtkGeodesicMeshDeformation.cxx | 274 +++++++ .../src/bbitkvtkGeodesicMeshDeformation.h | 75 ++ ...astMarchingQuadEdgeMeshFilterResultsBase.h | 193 +++++ ...tMarchingQuadEdgeMeshFilterResultsBase.hxx | 700 ++++++++++++++++++ 4 files changed, 1242 insertions(+) create mode 100644 packages/itkvtk/src/bbitkvtkGeodesicMeshDeformation.cxx create mode 100644 packages/itkvtk/src/bbitkvtkGeodesicMeshDeformation.h create mode 100644 packages/itkvtk/src/itkFastMarchingQuadEdgeMeshFilterResultsBase.h create mode 100644 packages/itkvtk/src/itkFastMarchingQuadEdgeMeshFilterResultsBase.hxx diff --git a/packages/itkvtk/src/bbitkvtkGeodesicMeshDeformation.cxx b/packages/itkvtk/src/bbitkvtkGeodesicMeshDeformation.cxx new file mode 100644 index 0000000..287c0e6 --- /dev/null +++ b/packages/itkvtk/src/bbitkvtkGeodesicMeshDeformation.cxx @@ -0,0 +1,274 @@ +//===== +// 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 +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 lstCenter = bbGetInputCenter(); + double s = bbGetInputS(); + bool ok = false; + using MeshType = itk::QuadEdgeMesh; + + //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; + 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; + 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; + 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 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 + + diff --git a/packages/itkvtk/src/bbitkvtkGeodesicMeshDeformation.h b/packages/itkvtk/src/bbitkvtkGeodesicMeshDeformation.h new file mode 100644 index 0000000..1ba0448 --- /dev/null +++ b/packages/itkvtk/src/bbitkvtkGeodesicMeshDeformation.h @@ -0,0 +1,75 @@ +//===== +// 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); + BBTK_DECLARE_INPUT(Direction, std::vector); + //BBTK_DECLARE_OUTPUT(Out,double); + BBTK_PROCESS(Process); + void Process(); + + using MeshType = itk::QuadEdgeMesh; + using FastMarchingType = itk::FastMarchingQuadEdgeMeshFilterResultsBase; + + long EdgeIdBack; + std::vector voiIdPoints; + std::vector 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,""); +BBTK_INPUT(GeodesicMeshDeformation,Direction,"(default [1,0,0]) [X,Y,Z]",std::vector,""); +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__ + diff --git a/packages/itkvtk/src/itkFastMarchingQuadEdgeMeshFilterResultsBase.h b/packages/itkvtk/src/itkFastMarchingQuadEdgeMeshFilterResultsBase.h new file mode 100644 index 0000000..a4d3b29 --- /dev/null +++ b/packages/itkvtk/src/itkFastMarchingQuadEdgeMeshFilterResultsBase.h @@ -0,0 +1,193 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#ifndef itkFastMarchingQuadEdgeMeshFilterBase_h +#define itkFastMarchingQuadEdgeMeshFilterBase_h + +#include "itkFastMarchingBase.h" +#include "itkFastMarchingTraits.h" + +namespace itk +{ +/** + *\class FastMarchingQuadEdgeMeshFilterBase + \brief Fast Marching Method on QuadEdgeMesh + + The speed function is specified by the input mesh. Data associated to each + point is considered as the speed function. The speed function is set using + the method SetInput(). + + If the speed function is constant and of value one, fast marching results is + an approximate geodesic function from the initial alive points. + + Implementation of this class is based on + "Fast Marching Methods on Triangulated Domains", Kimmel, R., and Sethian, J.A., + Proc. Nat. Acad. Sci., 95, pp. 8341-8435, 1998. + + \ingroup ITKFastMarching +*/ +template +//2023-04-27-PG class ITK_TEMPLATE_EXPORT FastMarchingQuadEdgeMeshFilterBase : public FastMarchingBase +class ITK_TEMPLATE_EXPORT FastMarchingQuadEdgeMeshFilterResultsBase : public FastMarchingBase +{ +public: + //ITK_DISALLOW_COPY_AND_MOVE(FastMarchingQuadEdgeMeshFilterBase); + ITK_DISALLOW_COPY_AND_MOVE(FastMarchingQuadEdgeMeshFilterResultsBase); + + //using Self = FastMarchingQuadEdgeMeshFilterBase; + using Self =FastMarchingQuadEdgeMeshFilterResultsBase; + using Superclass = FastMarchingBase; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + using Traits = typename Superclass::Traits; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + //itkTypeMacro(FastMarchingQuadEdgeMeshFilterBase, FastMarchingBase); + itkTypeMacro(FastMarchingQuadEdgeMeshFilterResultsBase, FastMarchingBase); + + using InputMeshType = typename Superclass::InputDomainType; + using InputMeshPointer = typename Superclass::InputDomainPointer; + using InputPixelType = typename Superclass::InputPixelType; + using InputPointType = typename InputMeshType::PointType; + using InputPointIdentifierType = typename InputMeshType::PointIdentifier; + + using OutputMeshType = typename Superclass::OutputDomainType; + using OutputMeshPointer = typename Superclass::OutputDomainPointer; + using OutputPixelType = typename Superclass::OutputPixelType; + using OutputPointType = typename OutputMeshType::PointType; + using OutputVectorType = typename OutputPointType::VectorType; + using OutputVectorRealType = typename OutputVectorType::RealValueType; + using OutputQEType = typename OutputMeshType::QEType; + using OutputPointIdentifierType = typename OutputMeshType::PointIdentifier; + using OutputPointsContainer = typename OutputMeshType::PointsContainer; + using OutputPointsContainerPointer = typename OutputPointsContainer::Pointer; + using OutputPointsContainerIterator = typename OutputPointsContainer::Iterator; + using OutputPointDataContainer = typename OutputMeshType::PointDataContainer; + using OutputPointDataContainerPointer = typename OutputPointDataContainer::Pointer; + + using OutputCellsContainer = typename OutputMeshType::CellsContainer; + using OutputCellsContainerPointer = typename OutputCellsContainer::Pointer; + using OutputCellsContainerConstIterator = typename OutputCellsContainer::ConstIterator; + using OutputCellType = typename OutputMeshType::CellType; + + + using NodeType = typename Traits::NodeType; + using NodePairType = typename Traits::NodePairType; + using NodePairContainerType = typename Traits::NodePairContainerType; + using NodePairContainerPointer = typename Traits::NodePairContainerPointer; + using NodePairContainerConstIterator = typename Traits::NodePairContainerConstIterator; + + // using NodeContainerType = typename Traits::NodeContainerType; + // using NodeContainerPointer = typename Traits::NodeContainerPointer; + // using NodeContainerConstIterator = typename Traits::NodeContainerConstIterator; + + using LabelType = typename Superclass::LabelType; + + using NodeLabelMapType = std::map; + using NodeLabelMapIterator = typename NodeLabelMapType::iterator; + using NodeLabelMapConstIterator = typename NodeLabelMapType::const_iterator; + +protected: + //FastMarchingQuadEdgeMeshFilterBase(); + FastMarchingQuadEdgeMeshFilterResultsBase(); + //~FastMarchingQuadEdgeMeshFilterBase() override = default; + ~FastMarchingQuadEdgeMeshFilterResultsBase() override = default; + + NodeLabelMapType m_Label; + + IdentifierType + GetTotalNumberOfNodes() const override; + + void + SetOutputValue(OutputMeshType * oMesh, const NodeType & iNode, const OutputPixelType & iValue) override; + + const OutputPixelType + GetOutputValue(OutputMeshType * oMesh, const NodeType & iNode) const override; + + unsigned char + GetLabelValueForGivenNode(const NodeType & iNode) const override; + + void + SetLabelValueForGivenNode(const NodeType & iNode, const LabelType & iLabel) override; + + void + UpdateNeighbors(OutputMeshType * oMesh, const NodeType & iNode) override; + + void + UpdateValue(OutputMeshType * oMesh, const NodeType & iNode) override; + + const OutputVectorRealType + Solve(OutputMeshType * oMesh, + const NodeType & iId, + const OutputPointType & iCurrentPoint, + const OutputVectorRealType & iF, + const NodeType & iId1, + const OutputPointType & iP1, + const bool & iIsFar1, + const OutputVectorRealType iVal1, + const NodeType & iId2, + const OutputPointType & iP2, + const bool & iIsFar2, + const OutputVectorRealType & iVal2) const; + + + const OutputVectorRealType + ComputeUpdate(const OutputVectorRealType & iVal1, + const OutputVectorRealType & iVal2, + const OutputVectorRealType & iNorm1, + const OutputVectorRealType & iSqNorm1, + const OutputVectorRealType & iNorm2, + const OutputVectorRealType & iSqNorm2, + const OutputVectorRealType & iDot, + const OutputVectorRealType & iF) const; + + bool + UnfoldTriangle(OutputMeshType * oMesh, + const OutputPointIdentifierType & iId, + const OutputPointType & iP, + const OutputPointIdentifierType & iId1, + const OutputPointType & iP1, + const OutputPointIdentifierType & iId2, + const OutputPointType & iP2, + OutputVectorRealType & oNorm, + OutputVectorRealType & oSqNorm, + OutputVectorRealType & oDot1, + OutputVectorRealType & oDot2, + OutputPointIdentifierType & oId) const; + + bool + CheckTopology(OutputMeshType * oMesh, const NodeType & iNode) override; + + void + InitializeOutput(OutputMeshType * oMesh) override; + +private: + const InputMeshType * m_InputMesh; + const InputMeshType * m_BackInputMesh; +}; +} // namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +//# include "itkFastMarchingQuadEdgeMeshFilterBase.hxx" +# include "itkFastMarchingQuadEdgeMeshFilterResultsBase.hxx" +#endif + +#endif // itkFastMarchingQuadEdgeMeshFilterBase_h diff --git a/packages/itkvtk/src/itkFastMarchingQuadEdgeMeshFilterResultsBase.hxx b/packages/itkvtk/src/itkFastMarchingQuadEdgeMeshFilterResultsBase.hxx new file mode 100644 index 0000000..5157413 --- /dev/null +++ b/packages/itkvtk/src/itkFastMarchingQuadEdgeMeshFilterResultsBase.hxx @@ -0,0 +1,700 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#ifndef itkFastMarchingQuadEdgeMeshFilterResultsBase_hxx +#define itkFastMarchingQuadEdgeMeshFilterResultsBase_hxx + +#include "itkFastMarchingQuadEdgeMeshFilterResultsBase.h" + +#include "itkMath.h" +#include "itkVector.h" +#include "itkMatrix.h" + +namespace itk +{ +template +FastMarchingQuadEdgeMeshFilterResultsBase::FastMarchingQuadEdgeMeshFilterResultsBase() + : Superclass() +{ + this->m_InputMesh = nullptr; + this->m_BackInputMesh = nullptr; +} + +template +IdentifierType +FastMarchingQuadEdgeMeshFilterResultsBase::GetTotalNumberOfNodes() const +{ + return this->GetInput()->GetNumberOfPoints(); +} + +template +void +FastMarchingQuadEdgeMeshFilterResultsBase::SetOutputValue(OutputMeshType * oMesh, + const NodeType & iNode, + const OutputPixelType & iValue) +{ + oMesh->SetPointData(iNode, iValue); +} + +template +const typename FastMarchingQuadEdgeMeshFilterResultsBase::OutputPixelType +FastMarchingQuadEdgeMeshFilterResultsBase::GetOutputValue(OutputMeshType * oMesh, + const NodeType & iNode) const +{ + OutputPixelType outputValue = NumericTraits::ZeroValue(); + oMesh->GetPointData(iNode, &outputValue); + return outputValue; +} + + +template +unsigned char +FastMarchingQuadEdgeMeshFilterResultsBase::GetLabelValueForGivenNode(const NodeType & iNode) const +{ + auto it = m_Label.find(iNode); + + if (it != m_Label.end()) + { + return it->second; + } + else + { + return Traits::Far; + } +} + +template +void +FastMarchingQuadEdgeMeshFilterResultsBase::SetLabelValueForGivenNode(const NodeType & iNode, + const LabelType & iLabel) +{ + m_Label[iNode] = iLabel; +} + +template +void +FastMarchingQuadEdgeMeshFilterResultsBase::UpdateNeighbors(OutputMeshType * oMesh, const NodeType & iNode) +{ + //printf("PG FastMarchingQuadEdgeMeshFilterResultsBase::UpdateNeighbors Start \n"); + OutputPointType p; + oMesh->GetPoint(iNode, &p); + + OutputQEType * qe = p.GetEdge(); + + if (qe) + { + OutputQEType * qe_it = qe; + this->m_InputMesh = this->GetInput(); + do + { + if (qe_it) + { + OutputPointIdentifierType neigh_id = qe_it->GetDestination(); + + const char label = this->GetLabelValueForGivenNode(neigh_id); + + if ((label != Traits::Alive) && (label != Traits::InitialTrial) && (label != Traits::Forbidden)) + { + this->UpdateValue(oMesh, neigh_id); + } + } + else + { + itkGenericExceptionMacro(<< "qe_it is nullptr"); + } + qe_it = qe_it->GetOnext(); + } while (qe_it != qe); + } + else + { + itkGenericExceptionMacro(<< "qe is nullptr"); + } + //printf("PG FastMarchingQuadEdgeMeshFilterResultsBase::UpdateNeighbors End \n"); +} + +template +void +FastMarchingQuadEdgeMeshFilterResultsBase::UpdateValue(OutputMeshType * oMesh, const NodeType & iNode) +{ + OutputPointType p; + oMesh->GetPoint(iNode, &p); + + InputPixelType F = NumericTraits::ZeroValue(); + //this->m_InputMesh->GetPointData(iNode, &F); + F = 1.; + if (F < 0.) + { + itkGenericExceptionMacro(<< "F < 0."); + } + + OutputQEType * qe = p.GetEdge(); + + OutputPixelType outputPixel = this->m_LargeValue; + + if (qe) + { + OutputQEType * qe_it = qe; + qe_it = qe_it->GetOnext(); + + do + { + OutputQEType * qe_it2 = qe_it->GetOnext(); + + if (qe_it2) + { + if (qe_it->GetLeft() != OutputMeshType::m_NoFace) + { + OutputPointIdentifierType id1 = qe_it->GetDestination(); + OutputPointIdentifierType id2 = qe_it2->GetDestination(); + + const auto label1 = static_cast(this->GetLabelValueForGivenNode(id1)); + const auto label2 = static_cast(this->GetLabelValueForGivenNode(id2)); + + bool IsFar1 = (label1 != Traits::Far); + bool IsFar2 = (label2 != Traits::Far); + + if (IsFar1 || IsFar2) + { + OutputPointType q1 = oMesh->GetPoint(id1); + OutputPointType q2 = oMesh->GetPoint(id2); + + auto val1 = static_cast(this->GetOutputValue(oMesh, id1)); + auto val2 = static_cast(this->GetOutputValue(oMesh, id2)); + + if (val1 > val2) + { + OutputPointType temp_pt(q1); + q1 = q2; + q2 = temp_pt; + + std::swap(val1, val2); + std::swap(IsFar1, IsFar2); + std::swap(id1, id2); + } + + const OutputVectorRealType temp = + this->Solve(oMesh, iNode, p, F, id1, q1, IsFar1, val1, id2, q2, IsFar2, val2); + + outputPixel = std::min(outputPixel, static_cast(temp)); + } + } + + qe_it = qe_it2; + } + else + { + // throw one exception here + itkGenericExceptionMacro(<< "qe_it2 is nullptr"); + } + } while (qe_it != qe); + + if (outputPixel < this->m_LargeValue) + { + this->SetOutputValue(oMesh, iNode, outputPixel); + + this->SetLabelValueForGivenNode(iNode, Traits::Trial); + + this->m_Heap.push(NodePairType(iNode, outputPixel)); + } + } + else + { + // throw one exception + itkGenericExceptionMacro(<< "qe_it is nullptr"); + } +} + +template +const typename FastMarchingQuadEdgeMeshFilterResultsBase::OutputVectorRealType +FastMarchingQuadEdgeMeshFilterResultsBase::Solve(OutputMeshType * oMesh, + const NodeType & iId, + const OutputPointType & iCurrentPoint, + const OutputVectorRealType & iF, + const NodeType & iId1, + const OutputPointType & iP1, + const bool & iIsFar1, + const OutputVectorRealType iVal1, + const NodeType & iId2, + const OutputPointType & iP2, + const bool & iIsFar2, + const OutputVectorRealType & iVal2) const +{ + OutputVectorType Edge1 = iP1 - iCurrentPoint; + OutputVectorType Edge2 = iP2 - iCurrentPoint; + + OutputVectorRealType sq_norm1 = Edge1.GetSquaredNorm(); + OutputVectorRealType norm1 = 0.; + + OutputVectorRealType epsilon = NumericTraits::epsilon(); + + if (sq_norm1 > epsilon) + { + norm1 = std::sqrt(sq_norm1); + + OutputVectorRealType inv_norm1 = 1. / norm1; + Edge1 *= inv_norm1; + } + + OutputVectorRealType sq_norm2 = Edge2.GetSquaredNorm(); + OutputVectorRealType norm2 = 0.; + + if (sq_norm2 > epsilon) + { + norm2 = std::sqrt(sq_norm2); + + OutputVectorRealType inv_norm2 = 1. / norm2; + Edge2 *= inv_norm2; + } + + if (!iIsFar1 && iIsFar2) + { + // only one point is a contributor + return iVal2 + norm2 * iF; + } + if (iIsFar1 && !iIsFar2) + { + // only one point is a contributor + return iVal1 + norm1 * iF; + } + + if (iIsFar1 && iIsFar2) + { + auto dot = static_cast(Edge1 * Edge2); + + if (dot >= 0.) + { + return ComputeUpdate(iVal1, iVal2, norm2, sq_norm2, norm1, sq_norm1, dot, iF); + } + else + { + OutputVectorRealType sq_norm3, norm3, dot1, dot2; + OutputPointIdentifierType new_id; + + bool unfolded = + UnfoldTriangle(oMesh, iId, iCurrentPoint, iId1, iP1, iId2, iP2, norm3, sq_norm3, dot1, dot2, new_id); + + if (unfolded) + { + OutputVectorRealType t_sq_norm3 = norm3 * norm3; + + auto val3 = static_cast(this->GetOutputValue(oMesh, new_id)); + OutputVectorRealType t1 = ComputeUpdate(iVal1, val3, norm3, t_sq_norm3, norm1, sq_norm1, dot1, iF); + OutputVectorRealType t2 = ComputeUpdate(iVal2, val3, norm3, t_sq_norm3, norm2, sq_norm2, dot2, iF); + + return std::min(t1, t2); + } + else + { + return ComputeUpdate(iVal1, iVal2, norm2, sq_norm2, norm1, sq_norm1, dot, iF); + } + } + } + + return this->m_LargeValue; +} + +template +const typename FastMarchingQuadEdgeMeshFilterResultsBase::OutputVectorRealType +FastMarchingQuadEdgeMeshFilterResultsBase::ComputeUpdate(const OutputVectorRealType & iVal1, + const OutputVectorRealType & iVal2, + const OutputVectorRealType & iNorm1, + const OutputVectorRealType & iSqNorm1, + const OutputVectorRealType & iNorm2, + const OutputVectorRealType & iSqNorm2, + const OutputVectorRealType & iDot, + const OutputVectorRealType & iF) const +{ + auto large_value = static_cast(this->m_LargeValue); + + OutputVectorRealType t = large_value; + + OutputVectorRealType CosAngle = iDot; + OutputVectorRealType SinAngle = std::sqrt(1. - iDot * iDot); + + OutputVectorRealType u = iVal2 - iVal1; + + OutputVectorRealType sq_u = u * u; + OutputVectorRealType f2 = iSqNorm1 + iSqNorm2 - 2. * iNorm1 * iNorm2 * CosAngle; + OutputVectorRealType f1 = iNorm2 * u * (iNorm1 * CosAngle - iNorm2); + OutputVectorRealType f0 = iSqNorm2 * (sq_u - iF * iF * iSqNorm1 * SinAngle * SinAngle); + + OutputVectorRealType delta = f1 * f1 - f0 * f2; + + OutputVectorRealType epsilon = NumericTraits::epsilon(); + + if (delta >= 0.) + { + if (itk::Math::abs(f2) > epsilon) + { + t = (-f1 - std::sqrt(delta)) / f2; + + // test if we must must choose the other solution + if ((t < u) || (iNorm2 * (t - u) / t < iNorm1 * CosAngle) || (iNorm1 / CosAngle < iNorm2 * (t - u) / t)) + { + t = (-f1 + std::sqrt(delta)) / f2; + } + } + else + { + if (itk::Math::abs(f1) > epsilon) + { + t = -f0 / f1; + } + else + { + t = -large_value; + } + } + } + else + { + t = -large_value; + } + + // choose the update from the 2 vertex only if upwind criterion is met + if ((u < t) && (iNorm1 * CosAngle < iNorm2 * (t - u) / t) && (iNorm2 * (t - u) / t < iNorm1 / CosAngle)) + { + return t + iVal1; + } + else + { + return std::min(iNorm2 * iF + iVal1, iNorm1 * iF + iVal2); + } +} + +template +bool +FastMarchingQuadEdgeMeshFilterResultsBase::UnfoldTriangle(OutputMeshType * oMesh, + const OutputPointIdentifierType & iId, + const OutputPointType & iP, + const OutputPointIdentifierType & iId1, + const OutputPointType & iP1, + const OutputPointIdentifierType & iId2, + const OutputPointType & iP2, + OutputVectorRealType & oNorm, + OutputVectorRealType & oSqNorm, + OutputVectorRealType & oDot1, + OutputVectorRealType & oDot2, + OutputPointIdentifierType & oId) const +{ + (void)iId; + + OutputVectorType Edge1 = iP1 - iP; + OutputVectorRealType Norm1 = Edge1.GetNorm(); + + OutputVectorRealType epsilon = NumericTraits::epsilon(); + + if (Norm1 > epsilon) + { + OutputVectorRealType inv_Norm = 1. / Norm1; + Edge1 *= inv_Norm; + } + + OutputVectorType Edge2 = iP2 - iP; + OutputVectorRealType Norm2 = Edge2.GetNorm(); + if (Norm2 > epsilon) + { + OutputVectorRealType inv_Norm = 1. / Norm2; + Edge2 *= inv_Norm; + } + + auto dot = static_cast(Edge1 * Edge2); + + // the equation of the lines defining the unfolding region + // [e.g. line 1 : {x ; =0} ] + using Vector2DType = Vector; + using Matrix2DType = Matrix; + + Vector2DType v1; + v1[0] = dot; + v1[1] = std::sqrt(1. - dot * dot); + + Vector2DType v2; + v2[0] = 1.; + v2[1] = 0.; + + Vector2DType x1; + x1[0] = Norm1; + x1[1] = 0.; + + Vector2DType x2 = Norm2 * v1; + + // keep track of the starting point + Vector2DType x_start1(x1); + Vector2DType x_start2(x2); + + OutputPointIdentifierType id1 = iId1; + OutputPointIdentifierType id2 = iId2; + + OutputQEType * qe = oMesh->FindEdge(id1, id2); + qe = qe->GetSym(); + + OutputPointType t_pt1 = iP1; + OutputPointType t_pt2 = iP2; + OutputPointType t_pt = t_pt1; + + unsigned int nNum = 0; + while (nNum < 50 && qe->GetLeft() != OutputMeshType::m_NoFace) + { + OutputQEType * qe_Lnext = qe->GetLnext(); + OutputPointIdentifierType t_id = qe_Lnext->GetDestination(); + + oMesh->GetPoint(t_id, &t_pt); + + Edge1 = t_pt2 - t_pt1; + Norm1 = Edge1.GetNorm(); + if (Norm1 > epsilon) + { + OutputVectorRealType inv_Norm = 1. / Norm1; + Edge1 *= inv_Norm; + } + + Edge2 = t_pt - t_pt1; + Norm2 = Edge2.GetNorm(); + if (Norm2 > epsilon) + { + OutputVectorRealType inv_Norm = 1. / Norm2; + Edge2 *= inv_Norm; + } + + /* compute the position of the new point x on the unfolding plane (via a rotation of -alpha on (x2-x1)/rNorm1 ) + | cos(alpha) sin(alpha)| + x = |-sin(alpha) cos(alpha)| * [x2-x1]*rNorm2/rNorm1 + x1 where cos(alpha)=dot + */ + Vector2DType vv; + if (Norm1 > epsilon) + { + vv = (x2 - x1) * Norm2 / Norm1; + } + else + { + vv.Fill(0.); + } + + dot = Edge1 * Edge2; + + Matrix2DType rotation; + rotation[0][0] = dot; + rotation[0][1] = std::sqrt(1. - dot * dot); + rotation[1][0] = -rotation[0][1]; + rotation[1][1] = dot; + + Vector2DType x = rotation * vv + x1; + + + /* compute the intersection points. + We look for x=x1+lambda*(x-x1) or x=x2+lambda*(x-x2) with =0, so */ + OutputVectorRealType lambda11 = -(x1 * v1) / ((x - x1) * v1); // left most + OutputVectorRealType lambda12 = -(x1 * v2) / ((x - x1) * v2); // right most + OutputVectorRealType lambda21 = -(x2 * v1) / ((x - x2) * v1); // left most + OutputVectorRealType lambda22 = -(x2 * v2) / ((x - x2) * v2); // right most + bool bIntersect11 = (lambda11 >= 0.) && (lambda11 <= 1.); + bool bIntersect12 = (lambda12 >= 0.) && (lambda12 <= 1.); + bool bIntersect21 = (lambda21 >= 0.) && (lambda21 <= 1.); + bool bIntersect22 = (lambda22 >= 0.) && (lambda22 <= 1.); + + if (bIntersect11 && bIntersect12) + { + qe = oMesh->FindEdge(id1, t_id); + qe = qe->GetSym(); + id2 = t_id; + t_pt2 = t_pt; + x2 = x; + } + else + { + if (bIntersect21 && bIntersect22) + { + qe = oMesh->FindEdge(id2, t_id); + qe = qe->GetSym(); + id1 = t_id; + t_pt1 = t_pt; + x1 = x; + } + else + { // bIntersect11 && !bIntersect12 && !bIntersect21 && bIntersect22 + + // that's it, we have found the point + oSqNorm = x.GetSquaredNorm(); + + if (oSqNorm > epsilon) + { + oNorm = std::sqrt(oSqNorm); + OutputVectorRealType temp_norm = x_start1.GetNorm(); + if (temp_norm > epsilon) + { + oDot1 = x * x_start1 / (oNorm * temp_norm); + } + else + { + oDot1 = 0.; + } + temp_norm = x_start2.GetNorm(); + if (temp_norm > epsilon) + { + oDot2 = x * x_start2 / (oNorm * temp_norm); + } + else + { + oDot2 = 0.; + } + } + else + { + oNorm = 0.; + oDot1 = 0.; + oDot2 = 0.; + } + + oId = t_id; + return true; + } + } + ++nNum; + } + + return false; +} + +template +bool +FastMarchingQuadEdgeMeshFilterResultsBase::CheckTopology(OutputMeshType * oMesh, const NodeType & iNode) +{ + (void)oMesh; + (void)iNode; + + return true; +} + +template +void +FastMarchingQuadEdgeMeshFilterResultsBase::InitializeOutput(OutputMeshType * oMesh) +{ + printf("PG FastMarchingQuadEdgeMeshFilterResultsBase::InitializeOutput Start \n"); + if(this->GetInput() != m_InputMesh || this->GetOutput() == NULL){ + this->CopyInputMeshToOutputMeshGeometry(); + } + else{ + printf("PG FastMarchingQuadEdgeMeshFilterResultsBase::InitializeOutput FLAG 0.5 \n"); + oMesh = this->GetOutput(); + } + printf("PG FastMarchingQuadEdgeMeshFilterResultsBase::InitializeOutput FLAG 1 \n"); + // Check that the input mesh is made of triangles + //We already now it is made of triangles + /** + { + OutputCellsContainerPointer cells = oMesh->GetCells(); + OutputCellsContainerConstIterator c_it = cells->Begin(); + OutputCellsContainerConstIterator c_end = cells->End(); + + while (c_it != c_end) + { + OutputCellType * cell = c_it.Value(); + + if (cell->GetNumberOfPoints() != 3) + { + itkGenericExceptionMacro(<< "Input mesh has non triangular faces"); + } + ++c_it; + } + }*/ + + OutputPointsContainerPointer points = oMesh->GetPoints(); + + OutputPointDataContainerPointer pointdata = OutputPointDataContainer::New(); + pointdata->Reserve(points->Size()); + + OutputPointsContainerIterator p_it = points->Begin(); + OutputPointsContainerIterator p_end = points->End(); + + while (p_it != p_end) + { + pointdata->SetElement(p_it->Index(), this->m_LargeValue); + ++p_it; + } + + oMesh->SetPointData(pointdata); + + m_Label.clear(); + + if (this->m_AlivePoints) + { + NodePairContainerConstIterator pointsIter = this->m_AlivePoints->Begin(); + NodePairContainerConstIterator pointsEnd = this->m_AlivePoints->End(); + while (pointsIter != pointsEnd) + { + // get node from alive points container + NodeType idx = pointsIter->Value().GetNode(); + + if (points->IndexExists(idx)) + { + OutputPixelType outputPixel = pointsIter->Value().GetValue(); + + this->SetLabelValueForGivenNode(idx, Traits::Alive); + this->SetOutputValue(oMesh, idx, outputPixel); + } + + ++pointsIter; + } + } + if (this->m_ForbiddenPoints) + { + NodePairContainerConstIterator pointsIter = this->m_ForbiddenPoints->Begin(); + NodePairContainerConstIterator pointsEnd = this->m_ForbiddenPoints->End(); + + OutputPixelType zero = NumericTraits::ZeroValue(); + + while (pointsIter != pointsEnd) + { + NodeType idx = pointsIter->Value().GetNode(); + + if (points->IndexExists(idx)) + { + this->SetLabelValueForGivenNode(idx, Traits::Forbidden); + this->SetOutputValue(oMesh, idx, zero); + } + + ++pointsIter; + } + } + if (this->m_TrialPoints) + { + NodePairContainerConstIterator pointsIter = this->m_TrialPoints->Begin(); + NodePairContainerConstIterator pointsEnd = this->m_TrialPoints->End(); + + while (pointsIter != pointsEnd) + { + NodeType idx = pointsIter->Value().GetNode(); + + if (points->IndexExists(idx)) + { + OutputPixelType outputPixel = pointsIter->Value().GetValue(); + + this->SetLabelValueForGivenNode(idx, Traits::InitialTrial); + this->SetOutputValue(oMesh, idx, outputPixel); + + this->m_Heap.push(pointsIter->Value()); + } + + ++pointsIter; + } + }printf("PG FastMarchingQuadEdgeMeshFilterResultsBase::InitializeOutput ENDED \n"); +} +} // namespace itk + +#endif // itkFastMarchingQuadEdgeMeshFilterResultsBase_hxx -- 2.47.1