[Insight-users] Re: IsInside() method in itk::MeshSpatialObject

[Yan Yang]

Hi Julien,

Many thanks for your reply. Would you like to have a look at my requirement?

I'd like to build a spatial relationship between a 3D mesh (tetrahedron) and
a 3D image. Since the mesh was extracted from CT images using a software, I
would like to map the mesh to the orignial image data in my program, so that
I can access all the voxels inside each elements. I don't know to build a
common physical coordinate system for these two objects.

The following is my code and running result, the test point [50, 20, 0]
should not in the tetrhedron with
node no      coordinate
6              65.8246 14.8065 0
1              40.8246 39.8065 25
3              40.8246 14.8065 25
5              65.8246 14.8065 25
however the return bool value of IsInside is true.

----------------------My
code------------------------------------------------------------------------
---
typedef itk::DefaultDynamicMeshTraits< float, 3, 3 > MeshTrait;
typedef itk::Mesh< MeshPixelType, 3, MeshTrait>  MeshType;
typedef MeshType::CellTraits       CellTraits;
typedef itk::CellInterface< float, CellTraits >  CellInterfaceType;
typedef itk::TetrahedronCell< CellInterfaceType >  TetraCellType;
typedef MeshType::PointType       PointType;
typedef MeshType::CellType       CellType;
typedef CellType::CellAutoPointer      CellAutoPointer;


  //access points in cells in a mesh
  typedef MeshType::CellsContainer::Iterator CellIterator;
  CellIterator cellIterator = m_mesh->GetCells()->Begin();
  CellIterator end = m_mesh->GetCells()->End();

   // Create a mesh with one tetrahedtron cell
  // in order to extract a bounding box for each cell
  MeshType::Pointer cellMesh = MeshType::New();
  cellMesh->SetCellsAllocationMethod(
  MeshType::CellsAllocatedDynamicallyCellByCell );
  CellAutoPointer oneCell;
  oneCell.TakeOwnership( new TetraCellType );
  CellAutoPointer testCell;
  testCell.TakeOwnership( new TetraCellType );

  typedef itk::MeshSpatialObject< MeshType > MeshSpatialObjectType;
  MeshSpatialObjectType::Pointer    cellMeshSpatialObject;

  unsigned int numberOfCells = m_mesh->GetNumberOfCells();

  for(unsigned int cellId = 0; cellId < numberOfCells; cellId++)
  {
 std::cout << std::endl << "Cell ID = " << cellId << std::endl;

 MeshType::CellType *cellptr1 = cellIterator.Value();
 TetraCellType *cellTetra1 = dynamic_cast<TetraCellType *>( cellptr1 );
 TetraCellType::PointIdIterator pit = cellTetra1->PointIdsBegin();
 TetraCellType::PointIdIterator pitEnd = cellTetra1->PointIdsEnd();


 int id = 0;
 while(pit != pitEnd)
 {
   MeshType::PointType pp;
   m_mesh->GetPoint( *pit-1, &pp );
   cout << *pit << "  ";
   cout << pp[0] << " " << pp[1] << " " << pp[2] << endl;


   //create a mesh with one tetrahedron cell
   cellMesh->SetPoint( *pit-1, pp );
   oneCell->SetPointId( id, *pit);

   pit++;
   id++;
 }

 cellMesh->SetCell( 0, oneCell );

 // print the point info in the mesh with one cell - tetrahedron
 cellMesh->GetCell( 0, testCell );
 TetraCellType::PointIdIterator pit2 = testCell->PointIdsBegin();
 TetraCellType::PointIdIterator pitEnd2 = testCell->PointIdsEnd();

 cout << "cell mesh------" << endl;
 while(pit2 != pitEnd2)
 {
   MeshType::PointType pp1;
   cellMesh->GetPoint( *pit2-1, &pp1 );

   cout << *pit2 << "  ";
   cout << pp1[0] << " " << pp1[1] << " " << pp1[2] << endl;

   pit2++;
 }

 // extract a bounding box for the mesh object - one tetra cell
 cellMeshSpatialObject = MeshSpatialObjectType::New();
 cellMeshSpatialObject->SetMesh( cellMesh );

 std::cout << "Cell mesh bounds " <<
  cellMeshSpatialObject->GetBoundingBox()->GetBounds() << std::endl;

 //cellMeshSpatialObject->Print(std::cout);


 MeshSpatialObjectType::PointType myPhysicalPoint;
 myPhysicalPoint[0] = 50;
 myPhysicalPoint[1] = 20;
 myPhysicalPoint[2] = 0;

 cout << "Is my physical point inside? : " <<
   cellMeshSpatialObject->IsInside( myPhysicalPoint ) << endl;
  }
----------------------------------------------------------------------------
-----

--------------------------running
result-----------------------------------------
Cell ID = 0
7  65.8246 39.8065 25
1  40.8246 39.8065 25
6  65.8246 14.8065 0
5  65.8246 14.8065 25
cell mesh------
7  65.8246 39.8065 25
1  40.8246 39.8065 25
6  65.8246 14.8065 0
5  65.8246 14.8065 25
Cell mesh bounds        [40.8246, 65.8246, 14.8065, 39.8065, 0, 25]
Is my physical point [50, 20, 0] inside? : 0

Cell ID = 1
6  65.8246 14.8065 0
3  40.8246 14.8065 25
1  40.8246 39.8065 25
2  40.8246 39.8065 0
cell mesh------
6  65.8246 14.8065 0
3  40.8246 14.8065 25
1  40.8246 39.8065 25
2  40.8246 39.8065 0
Cell mesh bounds        [40.8246, 65.8246, 14.8065, 39.8065, 0, 25]
Is my physical point [50, 20, 0] inside? : 0

Cell ID = 2
7  65.8246 39.8065 25
6  65.8246 14.8065 0
1  40.8246 39.8065 25
2  40.8246 39.8065 0
cell mesh------
7  65.8246 39.8065 25
6  65.8246 14.8065 0
1  40.8246 39.8065 25
2  40.8246 39.8065 0
Cell mesh bounds        [40.8246, 65.8246, 14.8065, 39.8065, 0, 25]
Is my physical point [50, 20, 0] inside? : 0

Cell ID = 3
6  65.8246 14.8065 0
2  40.8246 39.8065 0
4  40.8246 14.8065 0
3  40.8246 14.8065 25
cell mesh------
6  65.8246 14.8065 0
2  40.8246 39.8065 0
4  40.8246 14.8065 0
3  40.8246 14.8065 25
Cell mesh bounds        [40.8246, 65.8246, 14.8065, 39.8065, 0, 25]
Is my physical point [50, 20, 0] inside? : 0

Cell ID = 4
6  65.8246 14.8065 0
1  40.8246 39.8065 25
3  40.8246 14.8065 25
5  65.8246 14.8065 25
cell mesh------
6  65.8246 14.8065 0
1  40.8246 39.8065 25
3  40.8246 14.8065 25
5  65.8246 14.8065 25
Cell mesh bounds        [40.8246, 65.8246, 14.8065, 39.8065, 0, 25]
Is my physical point [50, 20, 0] inside? : 1

Cell ID = 5
6  65.8246 14.8065 0
2  40.8246 39.8065 0
7  65.8246 39.8065 25
8  65.8246 39.8065 0
cell mesh------
6  65.8246 14.8065 0
2  40.8246 39.8065 0
7  65.8246 39.8065 25
8  65.8246 39.8065 0
Cell mesh bounds        [40.8246, 65.8246, 14.8065, 39.8065, 0, 25]
Is my physical point [50, 20, 0] inside? : 0
----------------------------------------------------------------------------
---------------------

Thanks a lot.

Kind regards,

Yan


----- Original Message ----- 
From: "Julien Jomier" <jjomier at cs.unc.edu>
To: <y.yang at anglia.ac.uk>
Sent: Friday, December 09, 2005 6:29 PM
Subject: Re: IsInside() method in itk::MeshSpatialObject


> Hi Yan,
>
> Can you be more specific? and send me a sample code that shows the
problem?
> I cannot reproduce the problem you are describing.
>
> Let me know,
>
> Julien
>
> Yan Yang wrote:
> > Hi Julien,
> >
> > I am using ITK2.4.1. I am trying to query if a point is inside a
> > tetrahedron. I created a Mesh Spatial Object with only one tetrahedron
cell,
> > and used IsInside() method to query if a point within the cell's
bounding
> > box is inside the tetrahedron. But I am afraid there might be something
> > wrong with the return boolean value. I would be grateful if anyone could
> > tell me if the IsInside() function can check a point inside the mesh.
> >
> > Kind regards,
> >
> > Yan
> >
> >
> >
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20051210/486b29dd/attachment.html