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