]> Creatis software - bbtk.git/blob - packages/itkvtk/src/bbitkvtkGeodesicMeshDeformation.cxx
Clean Code
[bbtk.git] / packages / itkvtk / src / bbitkvtkGeodesicMeshDeformation.cxx
1 //===== 
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)
3 //===== 
4 #include "bbitkvtkGeodesicMeshDeformation.h"
5 #include "bbitkvtkPackage.h"
6 #include <vtkCellArrayIterator.h>
7 namespace bbitkvtk
8 {
9
10 BBTK_ADD_BLACK_BOX_TO_PACKAGE(itkvtk,GeodesicMeshDeformation)
11 BBTK_BLACK_BOX_IMPLEMENTATION(GeodesicMeshDeformation,bbtk::AtomicBlackBox);
12 //===== 
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)
14 //===== 
15 void GeodesicMeshDeformation::Process()
16 {
17 // THE MAIN PROCESSING METHOD BODY
18 //   Here we simply set the input 'In' value to the output 'Out'
19 //   And print out the output value
20 // INPUT/OUTPUT ACCESSORS ARE OF THE FORM :
21 //    void bbSet{Input|Output}NAME(const TYPE&)
22 //    const TYPE& bbGet{Input|Output}NAME() const 
23 //    Where :
24 //    * NAME is the name of the input/output
25 //      (the one provided in the attribute 'name' of the tag 'input')
26 //    * TYPE is the C++ type of the input/output
27 //      (the one provided in the attribute 'type' of the tag 'input')
28     std::vector<double> lstCenter = bbGetInputCenter();
29     double s   = bbGetInputS();
30     bool    ok  = false;
31     bool        pdChanged = false;
32     using MeshType = itk::QuadEdgeMesh<double, 3>;
33     
34     //Set up QuadEdge and filter every time polydata changes
35     if ((bbGetInputIn() != polydata) && (bbGetInputActive()==true) && (bbGetInputIn() != NULL))
36     {
37         pdChanged = true;
38         //Reset displacement
39         if(lstCenter.size() != NULL)
40         {
41                         backLstCenter[0] = lstCenter[0];
42                         backLstCenter[1] = lstCenter[1];
43                         backLstCenter[2] = lstCenter[2];
44         }
45         /**
46                 Create QuadEdge
47                 */
48                 using PointType = MeshType::PointType;
49                 //polydataBack
50                 polydata = bbGetInputIn();
51                 
52                 //points = points of input polydata
53                 vtkPoints* points=bbGetInputIn()->GetPoints();
54                 
55                 quadEdgeMesh = MeshType::New();
56                 double currentPoint[3] = {};
57                 PointType p;
58                 
59                 int numberOfPoints = points->GetNumberOfPoints();;
60                 for(int pointId = 0; pointId < numberOfPoints; pointId++){
61                         points->GetPoint(pointId, currentPoint);
62                         p[0] = currentPoint[0];
63                         p[1] = currentPoint[1];
64                         p[2] = currentPoint[2];
65                         quadEdgeMesh->SetPoint(pointId, p);
66                         //point data must be 1 on each point so that fast marching can calculate the distance.
67                         quadEdgeMesh->SetPointData(pointId, 1.);
68                 }
69
70                 //add cells to QuadEdge mesh
71                 auto cellIterator = vtk::TakeSmartPointer(bbGetInputIn()->GetPolys()->NewIterator());
72                 vtkIdList* cell;
73                 for(cellIterator->GoToFirstCell(); !cellIterator->IsDoneWithTraversal(); cellIterator->GoToNextCell()){
74                         cell = cellIterator->GetCurrentCell();
75                         quadEdgeMesh->AddFaceTriangle(cell->GetId(0), cell->GetId(1), cell->GetId(2));
76                 }
77                 
78                 using FastMarchingType = itk::FastMarchingQuadEdgeMeshFilterResultsBase<MeshType, MeshType>;    
79                 fmmFilter = FastMarchingType::New();
80                 fmmFilter->SetInput(quadEdgeMesh);
81                 
82                 //Filter and starting trial point removed because the range of point ids is not known when polydatas are split.
83                 
84                 //using NodePairType = FastMarchingType::NodePairType;
85                 //using NodePairContainerType = FastMarchingType::NodePairContainerType;
86
87                 //auto trial = NodePairContainerType::New();
88                 //NodePairType nodePair(0, 0.);
89                 //trial->push_back(nodePair);
90
91                 //fmmFilter->SetTrialPoints(trial);
92                 
93                 EdgeIdBack = bbGetInputEdgeId();
94
95                 using CriterionType = itk::FastMarchingThresholdStoppingCriterion<MeshType, MeshType>;
96                 auto criterion = CriterionType::New();
97                 sCurrent = 1;
98                 criterion->SetThreshold(1);
99                 fmmFilter->SetStoppingCriterion(criterion);
100                 fmmFilter->SetCollectPoints(true);
101                 fmmFilter->SetReleaseDataBeforeUpdateFlag(false);
102                 //fmmFilter->Update();
103
104                 //printf("PG GeodesicMeshDeformation::Process  Filter Created \n");
105     }
106     
107     if (bbGetInputTypeIn()==0) // direction
108     {
109         if (bbGetInputDirection().size()==3)
110         {
111             ok = !( (bbGetInputDirection()[0]==0) && (bbGetInputDirection()[1]==0) && (bbGetInputDirection()[2]==0) );
112         }
113     } // if TypeIn 0
114     
115     if (bbGetInputTypeIn()==1)  // center
116     {
117         ok = ( lstCenter.size()==3 );
118     } // if TypeIn 1
119
120     if ( (bbGetInputIn()!=NULL) && (ok==true) && (bbGetInputEdgeId()>=0) && (bbGetInputActive()==true) )
121     {
122         vtkPoints* points=bbGetInputIn()->GetPoints();
123         long    i,size=points->GetNumberOfPoints();
124         double  p[3];  // point
125         double  pb[3]; // point base
126         double  np[3]; // new point
127         double  sx,sy,sz;
128         sx = s*2;
129         sy = sx;
130         sz = sy;
131
132         double displcement_x = 0;
133         double displcement_y = 0;
134         double displcement_z = 0;
135         
136         if (bbGetInputTypeIn()==0) // Direction
137         {
138             displcement_x = bbGetInputDirection()[0];
139             displcement_y = bbGetInputDirection()[1];
140             displcement_z = bbGetInputDirection()[2];
141         } // if TypeIn 0 Direction
142                 
143                 printf(" EED GeodesicMeshDeformation::Process   %ld   %ld  -   %f %f %f \n", EdgeIdBack, bbGetInputEdgeId() , lstCenter[0],lstCenter[1],lstCenter[2] );
144                 
145         if (bbGetInputTypeIn()==1) // Center
146         {
147             if (EdgeIdBack==bbGetInputEdgeId() )
148             {
149                 displcement_x = (lstCenter[0]-backLstCenter[0])/1.0;
150                 displcement_y = (lstCenter[1]-backLstCenter[1])/1.0;
151                 displcement_z = (lstCenter[2]-backLstCenter[2])/1.0;
152             } // if EdgeIdBack!=bbGetInputEdgeId()
153             backLstCenter[0] = lstCenter[0];
154             backLstCenter[1] = lstCenter[1];
155             backLstCenter[2] = lstCenter[2];
156             
157         } // if TypeIn 1 Center
158         points->GetPoint( bbGetInputEdgeId() , pb );
159         if (EdgeIdBack!=bbGetInputEdgeId() || pdChanged)
160         {printf("PG GeodesicMeshDeformation::Process ENTERED \n");
161             EdgeIdBack = bbGetInputEdgeId();
162             
163                         backLstCenter[0] = lstCenter[0];
164                         backLstCenter[1] = lstCenter[1];
165                         backLstCenter[2] = lstCenter[2];
166                         
167                         using NodePairType = FastMarchingType::NodePairType;
168                         using NodePairContainerType = FastMarchingType::NodePairContainerType;
169                         
170                         auto processedPoints = NodePairContainerType::New();
171                         fmmFilter->SetProcessedPoints(processedPoints);
172
173                         auto trial = NodePairContainerType::New();
174
175                         //Seed point for fast marching with distance 0
176                         NodePairType nodePair(bbGetInputEdgeId(), 0.);
177                         trial->push_back(nodePair);
178                         
179                         fmmFilter->SetTrialPoints(trial);
180
181                         //threshold distance for fast marching method
182                         if(sCurrent != s)
183                         {                       
184                                 using CriterionType = itk::FastMarchingThresholdStoppingCriterion<MeshType, MeshType>;
185                                 auto criterion = CriterionType::New();
186                                 sCurrent = s;
187                                 criterion->SetThreshold(s*2);
188                                 fmmFilter->SetStoppingCriterion(criterion);
189                         }
190                         fmmFilter->Update();
191                         
192         } // if EdgeIdBack
193         
194         if ( !((displcement_x==0) &&(displcement_y==0) && (displcement_z==0)) )
195         {       
196                         if(fmmFilter != NULL)
197                         {
198                                 MeshType::PointDataContainer::Pointer pointData = fmmFilter->GetOutput()->GetPointData();
199                                 double pPD[3];
200                                 
201                                 FastMarchingType::NodePairContainerType::Iterator BegProcessedIt = fmmFilter->GetProcessedPoints()->Begin();
202                                 FastMarchingType::NodePairContainerType::Iterator EndProcessedIt = fmmFilter->GetProcessedPoints()->End();
203                                 
204                                 while (BegProcessedIt != EndProcessedIt)
205                                 {       
206                                         //point modification
207                                         points->GetPoint(BegProcessedIt.Value().GetNode(), pPD);
208                                         double weight = exp(-BegProcessedIt.Value().GetValue() * 1. / s);
209                                         if(weight > 0.0001){
210                                                 np[0] = pPD[0]+ weight*displcement_x;
211                                                 np[1] = pPD[1]+ weight*displcement_y;
212                                                 np[2] = pPD[2]+ weight*displcement_z;
213                                                 points->SetPoint(BegProcessedIt.Value().GetNode(), np);
214                                         }
215                                         ++BegProcessedIt;
216                                 }
217                                 points->Modified();
218                 bbGetInputIn()->Modified();
219                         }// if ffmFilter != NULL
220         } // if distance_x y z  != 0
221     } // In != NULL    ok    active
222 }
223 //===== 
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)
225 //===== 
226 void GeodesicMeshDeformation::bbUserSetDefaultValues()
227 {
228
229 //  SET HERE THE DEFAULT INPUT/OUTPUT VALUES OF YOUR BOX 
230 //    Here we initialize the input 'In' to 0
231         bbSetInputActive(true);
232         bbSetInputTypeIn(0);
233         bbSetInputIn(NULL);
234         std::vector<double> direction;
235         direction.push_back(1);
236         direction.push_back(0);
237         direction.push_back(0);
238         bbSetInputDirection(direction);
239         EdgeIdBack=-1;
240         bbSetInputEdgeId(EdgeIdBack);
241         bbSetInputS(10);
242
243         backLstCenter.push_back(0);
244         backLstCenter.push_back(0);
245         backLstCenter.push_back(0);
246         polydata = NULL;
247         quadEdgeMesh = NULL;
248         sCurrent = 0;
249         fmmFilter = NULL;
250         firstTime = true;
251 }
252 //===== 
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)
254 //===== 
255 void GeodesicMeshDeformation::bbUserInitializeProcessing()
256 {
257
258 //  THE INITIALIZATION METHOD BODY :
259 //    Here does nothing 
260 //    but this is where you should allocate the internal/output pointers 
261 //    if any 
262
263   
264 }
265 //===== 
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)
267 //===== 
268 void GeodesicMeshDeformation::bbUserFinalizeProcessing()
269 {
270
271 //  THE FINALIZATION METHOD BODY :
272 //    Here does nothing 
273 //    but this is where you should desallocate the internal/output pointers 
274 //    if any
275   
276 }
277 }
278 // EO namespace bbitkvtk
279
280