<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<HTML><HEAD>
<META http-equiv=Content-Type content="text/html; charset=iso-8859-1">
<META content="MSHTML 6.00.5730.11" name=GENERATOR>
<STYLE></STYLE>
</HEAD>
<BODY bgColor=#ffffff>
<DIV><FONT face=Arial size=2><FONT face="Times New Roman" size=3>Dear
Insight-Users,</FONT><BR></FONT></DIV>
<DIV><FONT face=Arial size=2>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.</FONT></DIV>
<DIV><FONT face=Arial size=2>Or the problem could be another?</FONT></DIV>
<DIV><FONT face=Arial size=2>Thank you in advance,</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2>Cecilia.</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2>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):</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2>#include "vtkPolyDataReader.h"<BR>#include
"vtkPolyDataWriter.h"<BR>#include "vtkAppendPolyData.h"<BR>#include
"vtkPolyData.h"<BR>#include "vtkOBBTree.h"<BR>#include
"vtkPolyDataConnectivityFilter.h"<BR>#include "vtkPoints.h"<BR>#include
"vtkFloatArray.h"<BR>#include "vtkPointData.h"<BR>#include "vtkDataSetWriter.h"
<BR>#include <vector><BR>#include <iostream><BR>#include
<stdio.h><BR>#include <math.h><BR>#include
<stdlib.h></FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2>bool isin(std::vector<int> v, int
a){<BR> for (int w = 0; w < v.size(); w++)<BR> if (v[w] == a)
return true;<BR> return false;<BR>}</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><BR><FONT face=Arial size=2>int main( int argc, char * argv[]
)<BR>{</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> float * point;<BR> FILE
*file;<BR> vtkPolyDataReader *CentersReader =
vtkPolyDataReader::New();<BR> vtkPolyDataReader *SurfaceReader =
vtkPolyDataReader::New();<BR> vtkPoints* outputPoints =
vtkPoints::New();<BR> vtkPolyData *append_centers =
vtkPolyData::New();<BR> vtkAppendPolyData *append_surface =
vtkAppendPolyData::New();<BR> vtkPolyDataConnectivityFilter *connectivity =
vtkPolyDataConnectivityFilter::New();<BR> vtkOBBTree *OBB =
vtkOBBTree::New();<BR> vtkPolyData *SavedData =
vtkPolyData::New();<BR> vtkFloatArray *FA =
vtkFloatArray::New();<BR> vtkPolyDataWriter *CentersWriter =
vtkPolyDataWriter::New();<BR> vtkPolyDataWriter *SurfaceWriter =
vtkPolyDataWriter::New();<BR> vtkDataSetWriter *DataSetWriter =
vtkDataSetWriter::New();<BR> <BR> if (argc !=
8){<BR> std::cout<<"Use: DeleteCenters.exe initialTimeStep
finalTimeStep inputCentersPath outputCentersPath inputSurfacePath
outputSurfacePath outputFilePath";<BR> return
0;<BR> }</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> char temp_a[200];<BR> char
temp_b[200];<BR> char temp_c[200];<BR> char temp_d[200];<BR> char
temp_file[200];</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> int initialTimeStep; //initial time step
<BR> ::sscanf(argv[1], "%d", &initialTimeStep);<BR> int
finalTimeStep; //final time step<BR> ::sscanf(argv[2], "%d",
&finalTimeStep);</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> printf
("**************************************");<BR> for (int m=initialTimeStep;
m<=finalTimeStep; m++){//cycle in time</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> printf("\n TIME =
%d\n",m);</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> if (m<10)
{<BR> sprintf(temp_a,"%s/Centers_t0%d.vtk",argv[3],
m);<BR> sprintf(temp_b,"%s/TotalIsosurface_t0%d.vtk",
argv[5],m);<BR> sprintf(temp_file,"%s/Informations_t0%d.txt",
argv[7],m);<BR> } else
{<BR> sprintf(temp_a,"%s/Centers_t%d.vtk",argv[3],
m);<BR> sprintf(temp_b,"%s/TotalIsosurface_t%d.vtk",
argv[5],m);<BR> sprintf(temp_file,"%s/Informations_t%d.txt",
argv[7],m);<BR> }</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial
size=2> file=fopen(temp_file,"w");<BR> if
(file==NULL){perror ("Error in file
creation.");exit(1);}<BR> fprintf(file, "\n TIME =
%d\n",m);</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial
size=2> CentersReader->SetFileName(temp_a);<BR> CentersReader->Update();</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial
size=2> SurfaceReader->SetFileName(temp_b);<BR> SurfaceReader->Update();</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> int
NumPoints=(CentersReader->GetOutput())->GetNumberOfPoints();<BR> fprintf(file,"\n
Number of centers found by the Hough Transform: %d\n",NumPoints);</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> <FONT
color=#ff0000>connectivity->SetInput(SurfaceReader->GetOutput());<BR> connectivity->ScalarConnectivityOn();</FONT></FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> std::vector<int>
AddedOrDouble;//lista dei centri e delle membrane già aggiunti o
doppi<BR> std::vector<int> Cancelled;//vettore dei centri
cancellati</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> int
added_centers_counter=0;</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> for (int i=0; i<NumPoints;
i++){//cycle of the connected surfaces (membranes)</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> if
(!isin(AddedOrDouble,i)){</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> float
range[2];<BR> range[0]=i;<BR> range[1]=i;</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> <FONT
color=#ff0000>connectivity->SetScalarRange(range);<BR> connectivity->Update();<BR> connectivity->Modified();</FONT></FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> //The class vtkOBBTree
determine whether a point is inside<BR> //or outside the
data used to build this OBB tree. The data <BR> //must be
a closed surface vtkPolyData data set. The return
<BR> //value is +1 if outside, -1 if inside, and 0 if
undecided. </FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> <FONT
color=#ff0000>OBB->SetDataSet(connectivity->GetOutput());<BR> OBB->Update();<BR> OBB->Modified();</FONT></FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial
size=2>/*DataSetWriter->SetInput(OBB->GetDataSet());<BR>DataSetWriter->SetFileName("C:/segm_t20_totale/modified/DataSet.vtk");<BR>DataSetWriter->Update();*/</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> int
result=0;<BR>// int added_centers_counter=0;</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> printf("\n Membrane %d,
",i);<BR> fprintf(file, "\n Membrane %d,
",i);</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> for (int k=0;
k<NumPoints; k++){//cycle of the centers found by the Hough
Transform</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> if
(!isin(AddedOrDouble,k)){</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial
size=2> point=(CentersReader->GetOutput())->GetPoints()->GetPoint(k);</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> if
(k==i){<BR> outputPoints->InsertPoint
(added_centers_counter,point[0], point[1],
point[2]);<BR> SavedData->DeepCopy(connectivity->GetOutput());<BR> int
nPoint=SavedData->GetNumberOfPoints();<BR> FA->SetNumberOfValues(nPoint);<BR> for
(float n=0;n<nPoint;n++
){<BR> FA->SetValue(n,added_centers_counter);<BR> }<BR> SavedData->GetPointData()->SetScalars(FA);<BR> append_surface->AddInput
(SavedData);<BR> added_centers_counter=added_centers_counter+1;<BR> AddedOrDouble.push_back(k);<BR> fprintf(file,"internal
centers: %d",k);<BR> }</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> <FONT
color=#ff0000>result=OBB->InsideOrOutside(point);</FONT></FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> if (k!=i
&&
result==-1){<BR> AddedOrDouble.push_back(k);<BR> fprintf(file,
",
%d",k);<BR> Cancelled.push_back(k);<BR> } <BR> } <BR> }<BR> fprintf(file,".\n");</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> if
(Cancelled.size()>0){<BR> for (int j=0;
j<Cancelled.size();
j++){<BR> fprintf(file, " Center %d,
cancelled.\n",
Cancelled[j]);<BR> }<BR> Cancelled.clear();<BR> }<BR> }<BR> }<BR> <BR> append_centers->SetPoints(outputPoints);<BR> fprintf(file,"\n
Number of centers after computation:
%d\n",append_centers->GetNumberOfPoints());<BR> fclose(file);</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> //Saving<BR> if (m<10)
{<BR> sprintf(temp_c,"%s/Centers_t0%d.vtk",
argv[4],m);<BR> sprintf(temp_d,"%s/TotalIsosurface_t0%d.vtk",
argv[6],m);<BR> } else
{<BR> sprintf(temp_c,"%s/Centers_t%d.vtk",
argv[4],m);<BR> sprintf(temp_d,"%s/TotalIsosurface_t%d.vtk",
argv[6],m);<BR> }</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial
size=2> CentersWriter->SetInput(append_centers);<BR> CentersWriter->SetFileName(temp_c);<BR> CentersWriter->Update();</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial
size=2> SurfaceWriter->SetInput(append_surface->GetOutput());<BR> SurfaceWriter->SetFileName(temp_d);<BR> SurfaceWriter->Update();</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> printf
("**************************************");<BR> <BR> }</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> return 0;</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> if (connectivity)
connectivity->Delete();<BR> if (outputPoints)
outputPoints->Delete();<BR> if (append_centers)
append_centers->Delete();<BR> if (append_surface)
append_surface->Delete();<BR> if (CentersReader)
CentersReader->Delete();<BR> if (SurfaceReader)
SurfaceReader->Delete();<BR> if (CentersWriter)
CentersWriter->Delete();<BR> if (SurfaceWriter)
SurfaceWriter->Delete();<BR> if (SavedData)
SavedData->Delete();<BR> if (FA)
FA->Delete();<BR> if (OBB) OBB->Delete();</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> }</FONT></DIV></BODY></HTML>