[Insight-users] How to know if a polydata is a closed surface

[Sylvain Jaume]

Hi Cecilia,

You can use vtkFeatureEdges to get the boundary edges (set feature edges
to off). Your polydata is closed if you get no boundary edges.

Sylvain

Cecilia Zanella wrote:
> Dear Insight-Users,
> I'm using the InsideOrOutside method of the vtkOBBTree class to know if 
> a point is inside a surface (polydata). It seems not to work in the 
> right way because sometime it returns -1 (internal point) even if the 
> point is outside. For this reason, I have the doubt to use a not closed 
> surface. Do anybody knows if there is a function to verify if the 
> surface I'm using is a closed polydata? It is an isosurface resulting 
> from the application of the vtkContourFilter and after extracted using 
> the vtkPolyDataConnectivityFilter.
> Or the problem could be another?
> Thank you in advance,
>  
> Cecilia.
>  
> Here is the code I implemented (the part in red is were I extract the 
> surface and use the vtkOBBTree class to verify if the point is inside 
> the surface itself):
>  
> #include "vtkPolyDataReader.h"
> #include "vtkPolyDataWriter.h"
> #include "vtkAppendPolyData.h"
> #include "vtkPolyData.h"
> #include "vtkOBBTree.h"
> #include "vtkPolyDataConnectivityFilter.h"
> #include "vtkPoints.h"
> #include "vtkFloatArray.h"
> #include "vtkPointData.h"
> #include "vtkDataSetWriter.h"
> #include <vector>
> #include <iostream>
> #include <stdio.h>
> #include <math.h>
> #include <stdlib.h>
>  
> bool isin(std::vector<int> v, int a){
>  for (int w = 0; w < v.size(); w++)
>   if (v[w] == a) return true;
>   return false;
> }
>  
> 
> int main( int argc, char * argv[] )
> {
>  
>  float * point;
>  FILE *file;
>  vtkPolyDataReader *CentersReader = vtkPolyDataReader::New();
>  vtkPolyDataReader *SurfaceReader = vtkPolyDataReader::New();
>  vtkPoints* outputPoints = vtkPoints::New();
>  vtkPolyData *append_centers = vtkPolyData::New();
>  vtkAppendPolyData *append_surface = vtkAppendPolyData::New();
>  vtkPolyDataConnectivityFilter *connectivity = 
> vtkPolyDataConnectivityFilter::New();
>  vtkOBBTree *OBB = vtkOBBTree::New();
>  vtkPolyData *SavedData = vtkPolyData::New();
>  vtkFloatArray *FA = vtkFloatArray::New();
>  vtkPolyDataWriter *CentersWriter = vtkPolyDataWriter::New();
>  vtkPolyDataWriter *SurfaceWriter = vtkPolyDataWriter::New();
>  vtkDataSetWriter *DataSetWriter = vtkDataSetWriter::New();
>  
>  if (argc != 8){
>   std::cout<<"Use: DeleteCenters.exe initialTimeStep finalTimeStep 
> inputCentersPath outputCentersPath inputSurfacePath outputSurfacePath 
> outputFilePath";
>   return 0;
>  }
>  
>  char temp_a[200];
>  char temp_b[200];
>  char temp_c[200];
>  char temp_d[200];
>  char temp_file[200];
>  
>  int initialTimeStep; //initial time step
>  ::sscanf(argv[1], "%d", &initialTimeStep);
>  int finalTimeStep; //final time step
>  ::sscanf(argv[2], "%d", &finalTimeStep);
>  
>  printf ("**************************************");
>  for (int m=initialTimeStep; m<=finalTimeStep; m++){//cycle in time
>  
>   printf("\n TIME = %d\n",m);
>  
>   if (m<10) {
>    sprintf(temp_a,"%s/Centers_t0%d.vtk",argv[3], m);
>    sprintf(temp_b,"%s/TotalIsosurface_t0%d.vtk", argv[5],m);
>    sprintf(temp_file,"%s/Informations_t0%d.txt", argv[7],m);
>   } else {
>    sprintf(temp_a,"%s/Centers_t%d.vtk",argv[3], m);
>    sprintf(temp_b,"%s/TotalIsosurface_t%d.vtk", argv[5],m);
>    sprintf(temp_file,"%s/Informations_t%d.txt", argv[7],m);
>   }
>  
>   file=fopen(temp_file,"w");
>   if (file==NULL){perror ("Error in file creation.");exit(1);}
>   fprintf(file, "\n TIME = %d\n",m);
>  
>   CentersReader->SetFileName(temp_a);
>   CentersReader->Update();
>  
>   SurfaceReader->SetFileName(temp_b);
>   SurfaceReader->Update();
>  
>   int NumPoints=(CentersReader->GetOutput())->GetNumberOfPoints();
>   fprintf(file,"\n Number of centers found by the Hough Transform: 
> %d\n",NumPoints);
>  
>   connectivity->SetInput(SurfaceReader->GetOutput());
>   connectivity->ScalarConnectivityOn();
>  
>   std::vector<int> AddedOrDouble;//lista dei centri e delle membrane già 
> aggiunti o doppi
>   std::vector<int> Cancelled;//vettore dei centri cancellati
>  
>   int added_centers_counter=0;
>  
>   for (int i=0; i<NumPoints; i++){//cycle of the connected surfaces 
> (membranes)
>  
>    if (!isin(AddedOrDouble,i)){
>  
>     float range[2];
>     range[0]=i;
>     range[1]=i;
>  
>     connectivity->SetScalarRange(range);
>     connectivity->Update();
>     connectivity->Modified();
>  
>     //The class vtkOBBTree determine whether a point is inside
>     //or outside the data used to build this OBB tree. The data
>     //must be a closed surface vtkPolyData data set. The return
>     //value is +1 if outside, -1 if inside, and 0 if undecided.
>  
>     OBB->SetDataSet(connectivity->GetOutput());
>     OBB->Update();
>     OBB->Modified();
>  
> /*DataSetWriter->SetInput(OBB->GetDataSet());
> DataSetWriter->SetFileName("C:/segm_t20_totale/modified/DataSet.vtk");
> DataSetWriter->Update();*/
>  
>     int result=0;
> //    int added_centers_counter=0;
>  
>     printf("\n Membrane %d, ",i);
>     fprintf(file, "\n Membrane %d, ",i);
>  
>     for (int k=0; k<NumPoints; k++){//cycle of the centers found by the 
> Hough Transform
>  
>      if (!isin(AddedOrDouble,k)){
>  
>       point=(CentersReader->GetOutput())->GetPoints()->GetPoint(k);
>  
>       if (k==i){
>        outputPoints->InsertPoint (added_centers_counter,point[0], 
> point[1], point[2]);
>        SavedData->DeepCopy(connectivity->GetOutput());
>        int nPoint=SavedData->GetNumberOfPoints();
>        FA->SetNumberOfValues(nPoint);
>        for (float n=0;n<nPoint;n++ ){
>         FA->SetValue(n,added_centers_counter);
>        }
>        SavedData->GetPointData()->SetScalars(FA);
>        append_surface->AddInput (SavedData);
>        added_centers_counter=added_centers_counter+1;
>        AddedOrDouble.push_back(k);
>        fprintf(file,"internal centers: %d",k);
>       }
>  
>       result=OBB->InsideOrOutside(point);
>  
>       if (k!=i && result==-1){
>        AddedOrDouble.push_back(k);
>        fprintf(file, ", %d",k);
>        Cancelled.push_back(k);
>       }    
>      }     
>     }
>     fprintf(file,".\n");
>  
>     if (Cancelled.size()>0){
>      for (int j=0; j<Cancelled.size(); j++){
>       fprintf(file, " Center %d, cancelled.\n", Cancelled[j]);
>      }
>      Cancelled.clear();
>     }
>    }
>   }
>   
>   append_centers->SetPoints(outputPoints);
>   fprintf(file,"\n Number of centers after computation: 
> %d\n",append_centers->GetNumberOfPoints());
>   fclose(file);
>  
>   //Saving
>   if (m<10) {
>    sprintf(temp_c,"%s/Centers_t0%d.vtk", argv[4],m);
>    sprintf(temp_d,"%s/TotalIsosurface_t0%d.vtk", argv[6],m);
>   } else {
>    sprintf(temp_c,"%s/Centers_t%d.vtk", argv[4],m);
>    sprintf(temp_d,"%s/TotalIsosurface_t%d.vtk", argv[6],m);
>   }
>  
>   CentersWriter->SetInput(append_centers);
>   CentersWriter->SetFileName(temp_c);
>   CentersWriter->Update();
>  
>   SurfaceWriter->SetInput(append_surface->GetOutput());
>   SurfaceWriter->SetFileName(temp_d);
>   SurfaceWriter->Update();
>  
>   printf ("**************************************");
>  
>   }
>  
>   return 0;
>  
>   if (connectivity) connectivity->Delete();
>   if (outputPoints) outputPoints->Delete();
>   if (append_centers) append_centers->Delete();
>   if (append_surface) append_surface->Delete();
>   if (CentersReader) CentersReader->Delete();
>   if (SurfaceReader) SurfaceReader->Delete();
>   if (CentersWriter) CentersWriter->Delete();
>   if (SurfaceWriter) SurfaceWriter->Delete();
>   if (SavedData) SavedData->Delete();
>   if (FA) FA->Delete();
>   if (OBB) OBB->Delete();
>  
>  }
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users