[Insight-users] How to know if a polydata is a closed surface
Mike Jackson
mike.jackson at imts.us
Mon Jan 29 09:10:18 EST 2007
IF your polydata is ALL triangles then you can use the following
formula:
numTriangles = (2 * numVertices) -4
This holds true for each enclosed set of triangles. In other words,
if you have say 5 different collections of enclosed triangulated
surfaces, then you need to run the formula for EACH enclosed surface.
--
Mike Jackson Senior Research Engineer
Innovative Management & Technology Services
On Jan 29, 2007, at 3:53 AM, 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20070129/81322b25/attachment.html
More information about the Insight-users
mailing list