[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