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 //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
  24 //    Where :
  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();
  31     bool    ok  = 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         //Reset displacement
  38         if(lstCenter.size() != NULL)
  39         {
  40                         backLstCenter[0] = lstCenter[0];
  41                         backLstCenter[1] = lstCenter[1];
  42                         backLstCenter[2] = lstCenter[2];
  43         }
  44         /**
  45                 Create QuadEdge
  46                 */
  47                 using PointType = MeshType::PointType;
  48                 //polydataBack
  49                 polydata = bbGetInputIn();
  50
  51                 //points = points of input polydata
  52                 vtkPoints* points=bbGetInputIn()->GetPoints();
  53
  54                 quadEdgeMesh = MeshType::New();
  55                 double currentPoint[3] = {};
  56                 PointType p;
  57
  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.);
  67                 }
  68
  69                 //add cells to QuadEdge mesh
  70                 auto cellIterator = vtk::TakeSmartPointer(bbGetInputIn()->GetPolys()->NewIterator());
  71                 vtkIdList* cell;
  72                 for(cellIterator->GoToFirstCell(); !cellIterator->IsDoneWithTraversal(); cellIterator->GoToNextCell()){
  73                         cell = cellIterator->GetCurrentCell();
  74                         quadEdgeMesh->AddFaceTriangle(cell->GetId(0), cell->GetId(1), cell->GetId(2));
  75                 }
  76
  77                 using FastMarchingType = itk::FastMarchingQuadEdgeMeshFilterResultsBase<MeshType, MeshType>;
  78                 fmmFilter = FastMarchingType::New();
  79                 fmmFilter->SetInput(quadEdgeMesh);
  80
  81                 //Filter and starting trial point removed because the range of point ids is not known when polydatas are split.
  82
  83                 //using NodePairType = FastMarchingType::NodePairType;
  84                 //using NodePairContainerType = FastMarchingType::NodePairContainerType;
  85
  86                 //auto trial = NodePairContainerType::New();
  87                 //NodePairType nodePair(0, 0.);
  88                 //trial->push_back(nodePair);
  89
  90                 //fmmFilter->SetTrialPoints(trial);
  91
  92                 EdgeIdBack = bbGetInputEdgeId();
  93
  94                 using CriterionType = itk::FastMarchingThresholdStoppingCriterion<MeshType, MeshType>;
  95                 auto criterion = CriterionType::New();
  96                 sCurrent = 1;
  97                 criterion->SetThreshold(1);
  98                 fmmFilter->SetStoppingCriterion(criterion);
  99                 fmmFilter->SetCollectPoints(true);
 100                 fmmFilter->SetReleaseDataBeforeUpdateFlag(false);
 101                 //fmmFilter->Update();
 102
 103                 //printf("PG GeodesicMeshDeformation::Process  Filter Created \n");
 104     }
 105
 106     if (bbGetInputTypeIn()==0) // direction
 107     {
 108         if (bbGetInputDirection().size()==3)
 109         {
 110             ok = !( (bbGetInputDirection()[0]==0) && (bbGetInputDirection()[1]==0) && (bbGetInputDirection()[2]==0) );
 111         }
 112     } // if TypeIn 0
 113
 114     if (bbGetInputTypeIn()==1)  // center
 115     {
 116         ok = ( lstCenter.size()==3 );
 117     } // if TypeIn 1
 118
 119     if ( (bbGetInputIn()!=NULL) && (ok==true) && (bbGetInputEdgeId()>=0) && (bbGetInputActive()==true) )
 120     {
 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
 126         double  sx,sy,sz;
 127         sx = s*2;
 128         sy = sx;
 129         sz = sy;
 130
 131         double displcement_x = 0;
 132         double displcement_y = 0;
 133         double displcement_z = 0;
 134
 135         if (bbGetInputTypeIn()==0) // Direction
 136         {
 137             displcement_x = bbGetInputDirection()[0];
 138             displcement_y = bbGetInputDirection()[1];
 139             displcement_z = bbGetInputDirection()[2];
 140         } // if TypeIn 0 Direction
 141
 142                 printf(" EED GeodesicMeshDeformation::Process   %ld   %ld  -   %f %f %f \n", EdgeIdBack, bbGetInputEdgeId() , lstCenter[0],lstCenter[1],lstCenter[2] );
 143
 144         if (bbGetInputTypeIn()==1) // Center
 145         {
 146             if (EdgeIdBack==bbGetInputEdgeId() )
 147             {
 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];
 155
 156         } // if TypeIn 1 Center
 157         points->GetPoint( bbGetInputEdgeId() , pb );
 158         if (EdgeIdBack!=bbGetInputEdgeId() )
 159         {
 160             EdgeIdBack = bbGetInputEdgeId();
 161
 162                         backLstCenter[0] = lstCenter[0];
 163                         backLstCenter[1] = lstCenter[1];
 164                         backLstCenter[2] = lstCenter[2];
 165
 166                         using NodePairType = FastMarchingType::NodePairType;
 167                         using NodePairContainerType = FastMarchingType::NodePairContainerType;
 168
 169                         auto processedPoints = NodePairContainerType::New();
 170                         fmmFilter->SetProcessedPoints(processedPoints);
 171
 172                         auto trial = NodePairContainerType::New();
 173
 174                         //Seed point for fast marching with distance 0
 175                         NodePairType nodePair(bbGetInputEdgeId(), 0.);
 176                         trial->push_back(nodePair);
 177
 178                         fmmFilter->SetTrialPoints(trial);
 179
 180                         //threshold distance for fast marching method
 181                         if(sCurrent != s)
 182                         {
 183                                 using CriterionType = itk::FastMarchingThresholdStoppingCriterion<MeshType, MeshType>;
 184                                 auto criterion = CriterionType::New();
 185                                 sCurrent = s;
 186                                 criterion->SetThreshold(s*2);
 187                                 fmmFilter->SetStoppingCriterion(criterion);
 188                         }
 189                         fmmFilter->Update();
 190
 191         } // if EdgeIdBack
 192
 193         if ( !((displcement_x==0) &&(displcement_y==0) && (displcement_z==0)) )
 194         {
 195                         if(fmmFilter != NULL)
 196                         {
 197                                 MeshType::PointDataContainer::Pointer pointData = fmmFilter->GetOutput()->GetPointData();
 198                                 double pPD[3];
 199
 200                                 FastMarchingType::NodePairContainerType::Iterator BegProcessedIt = fmmFilter->GetProcessedPoints()->Begin();
 201                                 FastMarchingType::NodePairContainerType::Iterator EndProcessedIt = fmmFilter->GetProcessedPoints()->End();
 202
 203                                 while (BegProcessedIt != EndProcessedIt)
 204                                 {
 205                                         //point modification
 206                                         points->GetPoint(BegProcessedIt.Value().GetNode(), pPD);
 207                                         double weight = exp(-BegProcessedIt.Value().GetValue() * 1. / s);
 208                                         if(weight > 0.0001){
 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);
 213                                         }
 214                                         ++BegProcessedIt;
 215                                 }
 216                                 points->Modified();
 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");
 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