#3501 Geodesic Deformation

[bbtk.git] / packages / itkvtk / src / bbitkvtkGeodesicMeshDeformation.cxx

diff --git a/packages/itkvtk/src/bbitkvtkGeodesicMeshDeformation.cxx b/packages/itkvtk/src/bbitkvtkGeodesicMeshDeformation.cxx
index 287c0e695f26474175e1a1fdf930831598f9033f..98559509fcff7fbd279be8b342c6a87be2b9c03e 100644 (file)
--- a/packages/itkvtk/src/bbitkvtkGeodesicMeshDeformation.cxx
+++ b/packages/itkvtk/src/bbitkvtkGeodesicMeshDeformation.cxx
@@ -14,7 +14,6 @@ BBTK_BLACK_BOX_IMPLEMENTATION(GeodesicMeshDeformation,bbtk::AtomicBlackBox);
 //===== 
 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
@@ -27,15 +26,21 @@ printf("PG GeodesicMeshDeformation::Process START \n");
 //    * 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>;
+    double  s           = bbGetInputS();
+    bool    ok          = false;
+    bool       pdChanged   = false;
+    using MeshType      = itk::QuadEdgeMesh<double, 3>;
+    std::vector<double> deformInfo;
+       bbSetOutputOut(deformInfo);
+       //std::vector<double> displacementVector;
     
     //Set up QuadEdge and filter every time polydata changes
     if ((bbGetInputIn() != polydata) && (bbGetInputActive()==true) && (bbGetInputIn() != NULL))
     {
+       pdChanged = true;
        //Reset displacement
-       if(lstCenter.size() != NULL){
+       if(lstCenter.size() != NULL)
+       {
                        backLstCenter[0] = lstCenter[0];
                        backLstCenter[1] = lstCenter[1];
                        backLstCenter[2] = lstCenter[2];
@@ -64,28 +69,31 @@ printf("PG GeodesicMeshDeformation::Process START \n");
                        //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;
+               //Filter and starting trial point removed because the range of point ids is not known when polydatas are split.
                
-               auto trial = NodePairContainerType::New();
-               NodePairType nodePair(0, 0.);
-               trial->push_back(nodePair);
+               //using NodePairType = FastMarchingType::NodePairType;
+               //using NodePairContainerType = FastMarchingType::NodePairContainerType;
+
+               //auto trial = NodePairContainerType::New();
+               //NodePairType nodePair(0, 0.);
+               //trial->push_back(nodePair);
+
+               //fmmFilter->SetTrialPoints(trial);
                
-               fmmFilter->SetTrialPoints(trial);
+               EdgeIdBack = bbGetInputEdgeId();
 
                using CriterionType = itk::FastMarchingThresholdStoppingCriterion<MeshType, MeshType>;
                auto criterion = CriterionType::New();
@@ -94,9 +102,9 @@ printf("PG GeodesicMeshDeformation::Process START \n");
                fmmFilter->SetStoppingCriterion(criterion);
                fmmFilter->SetCollectPoints(true);
                fmmFilter->SetReleaseDataBeforeUpdateFlag(false);
-               fmmFilter->Update();
+               //fmmFilter->Update();
 
-               printf("PG GeodesicMeshDeformation::Process  Filter Created \n");
+               //printf("PG GeodesicMeshDeformation::Process  Filter Created \n");
     }
     
     if (bbGetInputTypeIn()==0) // direction
@@ -119,10 +127,10 @@ printf("PG GeodesicMeshDeformation::Process START \n");
         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  sx,sy,sz;
+//        sx = s*2;
+//        sy = s*2;
+//        sz = s*2;
 
         double displcement_x = 0;
         double displcement_y = 0;
@@ -135,7 +143,7 @@ printf("PG GeodesicMeshDeformation::Process START \n");
             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] );
+               //printf(" EED GeodesicMeshDeformation::Process   %ld   %ld  -   %f %f %f \n", EdgeIdBack, bbGetInputEdgeId() , lstCenter[0],lstCenter[1],lstCenter[2] );
                
         if (bbGetInputTypeIn()==1) // Center
         {
@@ -151,8 +159,8 @@ printf("PG GeodesicMeshDeformation::Process START \n");
             
         } // if TypeIn 1 Center
         points->GetPoint( bbGetInputEdgeId() , pb );
-        if (EdgeIdBack!=bbGetInputEdgeId() )
-        {
+        if (EdgeIdBack!=bbGetInputEdgeId() || pdChanged)
+        {printf("PG GeodesicMeshDeformation::Process ENTERED \n");
             EdgeIdBack = bbGetInputEdgeId();
             
                        backLstCenter[0] = lstCenter[0];
@@ -174,11 +182,12 @@ printf("PG GeodesicMeshDeformation::Process START \n");
                        fmmFilter->SetTrialPoints(trial);
 
                        //threshold distance for fast marching method
-                       if(sCurrent != s){                      
+                       if(sCurrent != s)
+                       {                       
                                using CriterionType = itk::FastMarchingThresholdStoppingCriterion<MeshType, MeshType>;
                                auto criterion = CriterionType::New();
                                sCurrent = s;
-                               criterion->SetThreshold(s*2);
+                               criterion->SetThreshold(s*4);
                                fmmFilter->SetStoppingCriterion(criterion);
                        }
                        fmmFilter->Update();
@@ -187,7 +196,8 @@ printf("PG GeodesicMeshDeformation::Process START \n");
         
         if ( !((displcement_x==0) &&(displcement_y==0) && (displcement_z==0)) )
         {      
-                       if(fmmFilter != NULL){
+                       if(fmmFilter != NULL)
+                       {
                                MeshType::PointDataContainer::Pointer pointData = fmmFilter->GetOutput()->GetPointData();
                                double pPD[3];
                                
@@ -207,12 +217,16 @@ printf("PG GeodesicMeshDeformation::Process START \n");
                                        }
                                        ++BegProcessedIt;
                                }
+                               //double directionMoved[3] = {lstCenter[0]-displcement_x, lstCenter[1]-displcement_y, lstCenter[2]-displcement_z};
+                               //vtkMath::Normalize(directionMoved);
+                               std::vector<double> info{lstCenter[0],lstCenter[1],lstCenter[2], bbGetInputDirection()[0] , bbGetInputDirection()[1], bbGetInputDirection()[2],(double) bbGetInputEdgeId(), s};
+                               bbSetOutputOut(info);
+                               cout << "info updated" << endl;
                                points->Modified();
-               bbGetInputIn()->Modified();
-                       }
+                       bbGetInputIn()->Modified();
+                       }// if ffmFilter != NULL
         } // 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)
@@ -233,6 +247,8 @@ void GeodesicMeshDeformation::bbUserSetDefaultValues()
        EdgeIdBack=-1;
        bbSetInputEdgeId(EdgeIdBack);
        bbSetInputS(10);
+       std::vector<double> OutputVect;
+       bbSetOutputOut(OutputVect);
 
        backLstCenter.push_back(0);
        backLstCenter.push_back(0);