<HTML><BODY style="word-wrap: break-word; -khtml-nbsp-mode: space; -khtml-line-break: after-white-space; ">IF your polydata is ALL triangles then you can use the following formula:<DIV><BR class="khtml-block-placeholder"></DIV><DIV>numTriangles = (2 * numVertices) -4</DIV><DIV><BR class="khtml-block-placeholder"></DIV><DIV>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.</DIV><DIV><BR class="khtml-block-placeholder"></DIV><DIV><BR><DIV> <SPAN class="Apple-style-span" style="border-collapse: separate; border-spacing: 0px 0px; color: rgb(0, 0, 0); font-family: Bitstream Vera Sans; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; text-align: auto; -khtml-text-decorations-in-effect: none; text-indent: 0px; -apple-text-size-adjust: auto; text-transform: none; orphans: 2; white-space: normal; widows: 2; word-spacing: 0px; "><DIV style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">-- </DIV><DIV style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">Mike Jackson Senior Research Engineer</DIV><DIV style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">Innovative Management & Technology Services</DIV><BR class="Apple-interchange-newline"></SPAN> </DIV><BR><DIV><DIV>On Jan 29, 2007, at 3:53 AM, Cecilia Zanella wrote:</DIV><BR class="Apple-interchange-newline"><BLOCKQUOTE type="cite"> <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><DIV style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">_______________________________________________</DIV><DIV style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">Insight-users mailing list</DIV><DIV style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; "><A href="mailto:Insight-users@itk.org">Insight-users@itk.org</A></DIV><DIV style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; "><A href="http://www.itk.org/mailman/listinfo/insight-users">http://www.itk.org/mailman/listinfo/insight-users</A></DIV> </BLOCKQUOTE></DIV><BR></DIV></BODY></HTML>