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