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

[Cecilia Zanella]

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();

 }
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20070129/a75f3078/attachment-0001.htm