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