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)
40 backLstCenter[0] = lstCenter[0];
41 backLstCenter[1] = lstCenter[1];
42 backLstCenter[2] = lstCenter[2];
47 using PointType = MeshType::PointType;
49 polydata = bbGetInputIn();
51 //points = points of input polydata
52 vtkPoints* points=bbGetInputIn()->GetPoints();
54 quadEdgeMesh = MeshType::New();
55 double currentPoint[3] = {};
58 int numberOfPoints = points->GetNumberOfPoints();;
59 for(int pointId = 0; pointId < numberOfPoints; pointId++){
60 points->GetPoint(pointId, currentPoint);
61 p[0] = currentPoint[0];
62 p[1] = currentPoint[1];
63 p[2] = currentPoint[2];
64 quadEdgeMesh->SetPoint(pointId, p);
65 //point data must be 1 on each point so that fast marching can calculate the distance.
66 quadEdgeMesh->SetPointData(pointId, 1.);
69 //add cells to QuadEdge mesh
70 auto cellIterator = vtk::TakeSmartPointer(bbGetInputIn()->GetPolys()->NewIterator());
72 for(cellIterator->GoToFirstCell(); !cellIterator->IsDoneWithTraversal(); cellIterator->GoToNextCell()){
73 cell = cellIterator->GetCurrentCell();
74 quadEdgeMesh->AddFaceTriangle(cell->GetId(0), cell->GetId(1), cell->GetId(2));
77 using FastMarchingType = itk::FastMarchingQuadEdgeMeshFilterResultsBase<MeshType, MeshType>;
78 fmmFilter = FastMarchingType::New();
79 fmmFilter->SetInput(quadEdgeMesh);
81 //Filter and starting trial point removed because the range of point ids is not known when polydatas are split.
83 //using NodePairType = FastMarchingType::NodePairType;
84 //using NodePairContainerType = FastMarchingType::NodePairContainerType;
86 //auto trial = NodePairContainerType::New();
87 //NodePairType nodePair(0, 0.);
88 //trial->push_back(nodePair);
90 //fmmFilter->SetTrialPoints(trial);
92 EdgeIdBack = bbGetInputEdgeId();
94 using CriterionType = itk::FastMarchingThresholdStoppingCriterion<MeshType, MeshType>;
95 auto criterion = CriterionType::New();
97 criterion->SetThreshold(1);
98 fmmFilter->SetStoppingCriterion(criterion);
99 fmmFilter->SetCollectPoints(true);
100 fmmFilter->SetReleaseDataBeforeUpdateFlag(false);
101 //fmmFilter->Update();
103 //printf("PG GeodesicMeshDeformation::Process Filter Created \n");
106 if (bbGetInputTypeIn()==0) // direction
108 if (bbGetInputDirection().size()==3)
110 ok = !( (bbGetInputDirection()[0]==0) && (bbGetInputDirection()[1]==0) && (bbGetInputDirection()[2]==0) );
114 if (bbGetInputTypeIn()==1) // center
116 ok = ( lstCenter.size()==3 );
119 if ( (bbGetInputIn()!=NULL) && (ok==true) && (bbGetInputEdgeId()>=0) && (bbGetInputActive()==true) )
121 vtkPoints* points=bbGetInputIn()->GetPoints();
122 long i,size=points->GetNumberOfPoints();
123 double p[3]; // point
124 double pb[3]; // point base
125 double np[3]; // new point
131 double displcement_x = 0;
132 double displcement_y = 0;
133 double displcement_z = 0;
135 if (bbGetInputTypeIn()==0) // Direction
137 displcement_x = bbGetInputDirection()[0];
138 displcement_y = bbGetInputDirection()[1];
139 displcement_z = bbGetInputDirection()[2];
140 } // if TypeIn 0 Direction
142 printf(" EED GeodesicMeshDeformation::Process %ld %ld - %f %f %f \n", EdgeIdBack, bbGetInputEdgeId() , lstCenter[0],lstCenter[1],lstCenter[2] );
144 if (bbGetInputTypeIn()==1) // Center
146 if (EdgeIdBack==bbGetInputEdgeId() )
148 displcement_x = (lstCenter[0]-backLstCenter[0])/1.0;
149 displcement_y = (lstCenter[1]-backLstCenter[1])/1.0;
150 displcement_z = (lstCenter[2]-backLstCenter[2])/1.0;
151 } // if EdgeIdBack!=bbGetInputEdgeId()
152 backLstCenter[0] = lstCenter[0];
153 backLstCenter[1] = lstCenter[1];
154 backLstCenter[2] = lstCenter[2];
156 } // if TypeIn 1 Center
157 points->GetPoint( bbGetInputEdgeId() , pb );
158 if (EdgeIdBack!=bbGetInputEdgeId() )
160 EdgeIdBack = bbGetInputEdgeId();
162 backLstCenter[0] = lstCenter[0];
163 backLstCenter[1] = lstCenter[1];
164 backLstCenter[2] = lstCenter[2];
166 using NodePairType = FastMarchingType::NodePairType;
167 using NodePairContainerType = FastMarchingType::NodePairContainerType;
169 auto processedPoints = NodePairContainerType::New();
170 fmmFilter->SetProcessedPoints(processedPoints);
172 auto trial = NodePairContainerType::New();
174 //Seed point for fast marching with distance 0
175 NodePairType nodePair(bbGetInputEdgeId(), 0.);
176 trial->push_back(nodePair);
178 fmmFilter->SetTrialPoints(trial);
180 //threshold distance for fast marching method
183 using CriterionType = itk::FastMarchingThresholdStoppingCriterion<MeshType, MeshType>;
184 auto criterion = CriterionType::New();
186 criterion->SetThreshold(s*2);
187 fmmFilter->SetStoppingCriterion(criterion);
193 if ( !((displcement_x==0) &&(displcement_y==0) && (displcement_z==0)) )
195 if(fmmFilter != NULL)
197 MeshType::PointDataContainer::Pointer pointData = fmmFilter->GetOutput()->GetPointData();
200 FastMarchingType::NodePairContainerType::Iterator BegProcessedIt = fmmFilter->GetProcessedPoints()->Begin();
201 FastMarchingType::NodePairContainerType::Iterator EndProcessedIt = fmmFilter->GetProcessedPoints()->End();
203 while (BegProcessedIt != EndProcessedIt)
206 points->GetPoint(BegProcessedIt.Value().GetNode(), pPD);
207 double weight = exp(-BegProcessedIt.Value().GetValue() * 1. / s);
209 np[0] = pPD[0]+ weight*displcement_x;
210 np[1] = pPD[1]+ weight*displcement_y;
211 np[2] = pPD[2]+ weight*displcement_z;
212 points->SetPoint(BegProcessedIt.Value().GetNode(), np);
217 bbGetInputIn()->Modified();
218 }// if ffmFilter != NULL
219 } // if distance_x y z != 0
220 } // In != NULL ok active
221 // printf("PG GeodesicMeshDeformation::Process END \n");
224 // 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)
226 void GeodesicMeshDeformation::bbUserSetDefaultValues()
229 // SET HERE THE DEFAULT INPUT/OUTPUT VALUES OF YOUR BOX
230 // Here we initialize the input 'In' to 0
231 bbSetInputActive(true);
234 std::vector<double> direction;
235 direction.push_back(1);
236 direction.push_back(0);
237 direction.push_back(0);
238 bbSetInputDirection(direction);
240 bbSetInputEdgeId(EdgeIdBack);
243 backLstCenter.push_back(0);
244 backLstCenter.push_back(0);
245 backLstCenter.push_back(0);
253 // 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)
255 void GeodesicMeshDeformation::bbUserInitializeProcessing()
258 // THE INITIALIZATION METHOD BODY :
260 // but this is where you should allocate the internal/output pointers
266 // 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)
268 void GeodesicMeshDeformation::bbUserFinalizeProcessing()
271 // THE FINALIZATION METHOD BODY :
273 // but this is where you should desallocate the internal/output pointers
278 // EO namespace bbitkvtk