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

Julien Jomier jjomier at cs.unc.edu
Sun Dec 11 11:31:45 EST 2005


Hi Yan,

I've tested your code and I cannot reproduce the problem you are 
reporting. The test point (50,20,0) is reported as outside the mesh by 
the IsInside() function.

What version of ITK are you using?

Here a sample code I used. Can you test it and tell us if it's working 
or not?

--- Code starts here ---

typedef itk::DefaultDynamicMeshTraits< float , 3, 3 > MeshTrait;
typedef itk::Mesh<float,3,MeshTrait>                  MeshType;
typedef MeshType::CellTraits                          CellTraits;
typedef itk::CellInterface< float, CellTraits >       CellInterfaceType;
typedef itk::TetrahedronCell<CellInterfaceType>       TetraCellType;
typedef itk::MeshSpatialObject<MeshType>           MeshSpatialObjectType;
typedef MeshType::PointType                           PointType;
typedef MeshType::CellType                            CellType;
typedef CellType::CellAutoPointer                     CellAutoPointer;

// Create an itkMesh
MeshType::Pointer mesh = MeshType::New();

MeshType::CoordRepType testPointCoords[4][3]
     = { {65.8246,14.8065,0}, {40.8246,39.8065,25}, 
{40.8246,14.8065,25}, {65.8246,14.8065,25} };

unsigned long tetraPoints[4] = {0,1,2,3};

int i;
for(i=0; i < 4 ; ++i)
   {
   mesh->SetPoint(i, PointType(testPointCoords[i]));
   }

 
mesh->SetCellsAllocationMethod(MeshType::CellsAllocatedDynamicallyCellByCell 
);
CellAutoPointer testCell1;
testCell1.TakeOwnership(  new TetraCellType );
testCell1->SetPointIds(tetraPoints);
mesh->SetCell(0, testCell1 );

// Create the mesh Spatial Object
MeshSpatialObjectType::Pointer meshSO = MeshSpatialObjectType::New();
meshSO->SetMesh(mesh);

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

std::cout << outside << " is inside ?" << meshSO->IsInside(outside) << 
std::endl;

--- code ends here ---

best regards,

Julien

Yan Yang wrote:
> 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
>>>
>>>
>>>
> 
> 
> 
> 


More information about the Insight-users mailing list