2 // 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)
4 #include "bbitkvtkGeodesicMeshDeformation.h"
5 #include "bbitkvtkPackage.h"
6 #include <vtkCellArrayIterator.h>
10 BBTK_ADD_BLACK_BOX_TO_PACKAGE(itkvtk,GeodesicMeshDeformation)
11 BBTK_BLACK_BOX_IMPLEMENTATION(GeodesicMeshDeformation,bbtk::AtomicBlackBox);
13 // 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)
15 void GeodesicMeshDeformation::Process()
17 printf("PG GeodesicMeshDeformation::Process START \n");
18 // THE MAIN PROCESSING METHOD BODY
19 // Here we simply set the input 'In' value to the output 'Out'
20 // And print out the output value
21 // INPUT/OUTPUT ACCESSORS ARE OF THE FORM :
22 // void bbSet{Input|Output}NAME(const TYPE&)
23 // const TYPE& bbGet{Input|Output}NAME() const
25 // * NAME is the name of the input/output
26 // (the one provided in the attribute 'name' of the tag 'input')
27 // * TYPE is the C++ type of the input/output
28 // (the one provided in the attribute 'type' of the tag 'input')
29 std::vector<double> lstCenter = bbGetInputCenter();
30 double s = bbGetInputS();
32 using MeshType = itk::QuadEdgeMesh<double, 3>;
34 //Set up QuadEdge and filter every time polydata changes
35 if ((bbGetInputIn() != polydata) && (bbGetInputActive()==true) && (bbGetInputIn() != NULL))
38 if(lstCenter.size() != NULL){
39 backLstCenter[0] = lstCenter[0];
40 backLstCenter[1] = lstCenter[1];
41 backLstCenter[2] = lstCenter[2];
46 using PointType = MeshType::PointType;
48 polydata = bbGetInputIn();
50 //points = points of input polydata
51 vtkPoints* points=bbGetInputIn()->GetPoints();
53 quadEdgeMesh = MeshType::New();
54 double currentPoint[3] = {};
57 int numberOfPoints = points->GetNumberOfPoints();;
58 for(int pointId = 0; pointId < numberOfPoints; pointId++){
59 points->GetPoint(pointId, currentPoint);
60 p[0] = currentPoint[0];
61 p[1] = currentPoint[1];
62 p[2] = currentPoint[2];
63 quadEdgeMesh->SetPoint(pointId, p);
64 //point data must be 1 on each point so that fast marching can calculate the distance.
65 quadEdgeMesh->SetPointData(pointId, 1.);
67 printf("PG GeodesicMeshDeformation::Process Entered FLAG 4 \n");
68 //add cells to QuadEdge mesh
69 auto cellIterator = vtk::TakeSmartPointer(bbGetInputIn()->GetPolys()->NewIterator());
71 for(cellIterator->GoToFirstCell(); !cellIterator->IsDoneWithTraversal(); cellIterator->GoToNextCell()){
72 cell = cellIterator->GetCurrentCell();
73 quadEdgeMesh->AddFaceTriangle(cell->GetId(0), cell->GetId(1), cell->GetId(2));
75 printf("PG GeodesicMeshDeformation::Process FLAG 5\n");
77 using FastMarchingType = itk::FastMarchingQuadEdgeMeshFilterResultsBase<MeshType, MeshType>;
78 fmmFilter = FastMarchingType::New();
79 fmmFilter->SetInput(quadEdgeMesh);
81 using NodePairType = FastMarchingType::NodePairType;
82 using NodePairContainerType = FastMarchingType::NodePairContainerType;
84 auto trial = NodePairContainerType::New();
85 NodePairType nodePair(0, 0.);
86 trial->push_back(nodePair);
88 fmmFilter->SetTrialPoints(trial);
90 using CriterionType = itk::FastMarchingThresholdStoppingCriterion<MeshType, MeshType>;
91 auto criterion = CriterionType::New();
93 criterion->SetThreshold(1);
94 fmmFilter->SetStoppingCriterion(criterion);
95 fmmFilter->SetCollectPoints(true);
96 fmmFilter->SetReleaseDataBeforeUpdateFlag(false);
99 printf("PG GeodesicMeshDeformation::Process Filter Created \n");
102 if (bbGetInputTypeIn()==0) // direction
104 if (bbGetInputDirection().size()==3)
106 ok = !( (bbGetInputDirection()[0]==0) && (bbGetInputDirection()[1]==0) && (bbGetInputDirection()[2]==0) );
110 if (bbGetInputTypeIn()==1) // center
112 ok = ( lstCenter.size()==3 );
115 if ( (bbGetInputIn()!=NULL) && (ok==true) && (bbGetInputEdgeId()>=0) && (bbGetInputActive()==true) )
117 vtkPoints* points=bbGetInputIn()->GetPoints();
118 long i,size=points->GetNumberOfPoints();
119 double p[3]; // point
120 double pb[3]; // point base
121 double np[3]; // new point
127 double displcement_x = 0;
128 double displcement_y = 0;
129 double displcement_z = 0;
131 if (bbGetInputTypeIn()==0) // Direction
133 displcement_x = bbGetInputDirection()[0];
134 displcement_y = bbGetInputDirection()[1];
135 displcement_z = bbGetInputDirection()[2];
136 } // if TypeIn 0 Direction
138 printf(" EED GeodesicMeshDeformation::Process %ld %ld - %f %f %f \n", EdgeIdBack, bbGetInputEdgeId() , lstCenter[0],lstCenter[1],lstCenter[2] );
140 if (bbGetInputTypeIn()==1) // Center
142 if (EdgeIdBack==bbGetInputEdgeId() )
144 displcement_x = (lstCenter[0]-backLstCenter[0])/1.0;
145 displcement_y = (lstCenter[1]-backLstCenter[1])/1.0;
146 displcement_z = (lstCenter[2]-backLstCenter[2])/1.0;
147 } // if EdgeIdBack!=bbGetInputEdgeId()
148 backLstCenter[0] = lstCenter[0];
149 backLstCenter[1] = lstCenter[1];
150 backLstCenter[2] = lstCenter[2];
152 } // if TypeIn 1 Center
153 points->GetPoint( bbGetInputEdgeId() , pb );
154 if (EdgeIdBack!=bbGetInputEdgeId() )
156 EdgeIdBack = bbGetInputEdgeId();
158 backLstCenter[0] = lstCenter[0];
159 backLstCenter[1] = lstCenter[1];
160 backLstCenter[2] = lstCenter[2];
162 using NodePairType = FastMarchingType::NodePairType;
163 using NodePairContainerType = FastMarchingType::NodePairContainerType;
165 auto processedPoints = NodePairContainerType::New();
166 fmmFilter->SetProcessedPoints(processedPoints);
168 auto trial = NodePairContainerType::New();
170 //Seed point for fast marching with distance 0
171 NodePairType nodePair(bbGetInputEdgeId(), 0.);
172 trial->push_back(nodePair);
174 fmmFilter->SetTrialPoints(trial);
176 //threshold distance for fast marching method
178 using CriterionType = itk::FastMarchingThresholdStoppingCriterion<MeshType, MeshType>;
179 auto criterion = CriterionType::New();
181 criterion->SetThreshold(s*2);
182 fmmFilter->SetStoppingCriterion(criterion);
188 if ( !((displcement_x==0) &&(displcement_y==0) && (displcement_z==0)) )
190 if(fmmFilter != NULL){
191 MeshType::PointDataContainer::Pointer pointData = fmmFilter->GetOutput()->GetPointData();
194 FastMarchingType::NodePairContainerType::Iterator BegProcessedIt = fmmFilter->GetProcessedPoints()->Begin();
195 FastMarchingType::NodePairContainerType::Iterator EndProcessedIt = fmmFilter->GetProcessedPoints()->End();
197 while (BegProcessedIt != EndProcessedIt)
200 points->GetPoint(BegProcessedIt.Value().GetNode(), pPD);
201 double weight = exp(-BegProcessedIt.Value().GetValue() * 1. / s);
203 np[0] = pPD[0]+ weight*displcement_x;
204 np[1] = pPD[1]+ weight*displcement_y;
205 np[2] = pPD[2]+ weight*displcement_z;
206 points->SetPoint(BegProcessedIt.Value().GetNode(), np);
211 bbGetInputIn()->Modified();
213 } // if distance_x y z != 0
214 } // In != NULL ok active
215 printf("PG GeodesicMeshDeformation::Process END \n");
218 // 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)
220 void GeodesicMeshDeformation::bbUserSetDefaultValues()
223 // SET HERE THE DEFAULT INPUT/OUTPUT VALUES OF YOUR BOX
224 // Here we initialize the input 'In' to 0
225 bbSetInputActive(true);
228 std::vector<double> direction;
229 direction.push_back(1);
230 direction.push_back(0);
231 direction.push_back(0);
232 bbSetInputDirection(direction);
234 bbSetInputEdgeId(EdgeIdBack);
237 backLstCenter.push_back(0);
238 backLstCenter.push_back(0);
239 backLstCenter.push_back(0);
247 // 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)
249 void GeodesicMeshDeformation::bbUserInitializeProcessing()
252 // THE INITIALIZATION METHOD BODY :
254 // but this is where you should allocate the internal/output pointers
260 // 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)
262 void GeodesicMeshDeformation::bbUserFinalizeProcessing()
265 // THE FINALIZATION METHOD BODY :
267 // but this is where you should desallocate the internal/output pointers
272 // EO namespace bbitkvtk