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