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