[Insight-users] Fwd: Labelling artifacts in itk::TriangleMeshToBinaryImageFilter
Ramón Casero Cañas
rcasero at gmail.com
Fri Oct 11 12:52:28 EDT 2013
Dear all,
I'm using itk::TriangleMeshToBinaryImageFilter (ITK v4.3.1) to rasterize a
triangular mesh as a binary image.
The rasterization seems to produce many mis-labelled voxels inside the
mesh. I tried to attach two CSV files with the mesh vertices and triangles
to help reproduce the result, but it seems that it's too big for the
insight-users mailing list. I can provide them via email, though.
I attach a picture of the mesh as well as resulting slices 231 and 300, to
illustrate the problem.
The output image parameters I'm using are:
Resolution (spacing) = [2.539e-05, 2.539e-05, 2.441e-05]
Size = [415, 460, 900]
Origin = [0.006221, 0.008633, 0.000976]
Some things I have tried:
* Play with the Tolerance value. In this example the program sets it
to 2.4410e-06. I have tried reducing it to 1e-12, but the solution doesn't
improve.
* Scale the mesh vertex coordinates, and correspondingly the output image
voxel size and origin coordinates. The solution is the same.
I have programmed another function that doesn't suffer from these artifacts
(using the library CGAL, without ITK), but it's a lot slower than the ITK
solution, so it'd be great to make the ITK approach work
https://code.google.com/p/gerardus/source/browse/trunk/matlab/CgalToolbox/CgalInSurfaceTriangulation.cpp?r=1347
I'm running the ITK-based solution as a MEX function from Matlab. The full
code is here
https://code.google.com/p/gerardus/source/browse/trunk/matlab/ItkToolbox/ItkTriRasterization.cpp?r=1347
but the relevants part of the code to ITK are
/* ITK headers */
#include <itkImage.h>
#include <itkMesh.h>
#include <itkTriangleMeshToBinaryImageFilter.h>
#include <itkTriangleCell.h>
// type definitions
static const unsigned int Dimension = 3;
typedef uint8_T PixelType;
typedef float MeshDataType;
typedef itk::PointSet<MeshDataType, Dimension> PointSetType;
typedef PointSetType::CoordRepType CoordType;
typedef itk::Mesh<MeshDataType, Dimension> MeshType;
typedef itk::Image<PixelType, Dimension> ImageType;
typedef PointSetType::PointType PointType;
typedef MeshType::CellType CellType;
typedef itk::TriangleCell<CellType> TriangleType;
typedef CellType::CellAutoPointer CellAutoPointer;
typedef itk::TriangleMeshToBinaryImageFilter<MeshType, ImageType>
MeshFilterType;
/*
* mexFunction(): entry point for the mex function
*/
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
... // instantiate mesh
MeshType::Pointer mesh = MeshType::New();
// read vertices
... PointSetType::Pointer x = PointSetType::New(); ...
// populate mesh with vertices
mesh->SetPoints(x->GetPoints());
// read triangles
PointType triDef;
triDef.Fill(mxGetNaN());
for (mwIndex i = 0; i < nrowsTRI; ++i) {
PointType triangle = matlabImport->ReadRowVectorFromMatlab<CoordType,
PointType>(inTRI, i, triDef);
// create a triangle cell to read the vertex indices of the current
input triangle
CellAutoPointer cell;
cell.TakeOwnership(new TriangleType);
// assign to the 0, 1, 2 elements in the triangle cell the vertex
// indices that we have just read. Note that we have to substract
// 1 to convert Matlab's index convention 1, 2, 3, ... to C++
// convention 0, 1, 2, ...
cell->SetPointId(0, triangle[0] - 1);
cell->SetPointId(1, triangle[1] - 1);
cell->SetPointId(2, triangle[2] - 1);
// insert cell into the mesh
mesh->SetCell(i, cell);
}
... // get user input parameters for the output rasterization
... ImageType::SpacingType spacing = ...
... ImageType::SizeType size = ...
... ImageType::PointType origin = ... // instantiate rasterization
filter
MeshFilterType::Pointer meshFilter = MeshFilterType::New();
// smallest voxel side length ImageType::SpacingValueType minSpacing =
spacing[0]; for (mwIndex i = 1; i < Dimension; ++i) { minSpacing =
std::min(minSpacing, spacing[i]); }
// pass input parameters to the filter
meshFilter->SetInput(mesh);
meshFilter->SetSpacing(spacing);
meshFilter->SetSize(size);
meshFilter->SetOrigin(origin);
meshFilter->SetTolerance(minSpacing / 10.0);
meshFilter->SetInsideValue(1);
meshFilter->SetOutsideValue(0);
ImageType::IndexType start;
start.Fill(0);
meshFilter->SetIndex(start);
// convert image size from itk::Size format to std::vector<mwSize>
// so that we can use it in GraftItkImageOntoMatlab
std::vector<mwSize> sizeStdVector(Dimension);
for (unsigned int i = 0; i < Dimension; ++i) {
sizeStdVector[i] = size[i];
}
// graft ITK filter outputs onto Matlab outputs
matlabExport->GraftItkImageOntoMatlab<PixelType, Dimension>
(outIM, meshFilter->GetOutput(), sizeStdVector);
...
// run rasterization
meshFilter->Update();
...
}
Any ideas are welcome, thanks!
Best regards,
Ramon.
--
Dr. Ramón Casero Cañas
Oxford e-Research Centre (OeRC)
University of Oxford
7 Keble Rd
Oxford OX1 3QG
tlf +44 (0) 1865 610739
web http://www.cs.ox.ac.uk/people/Ramon.CaseroCanas
photos http://www.flickr.com/photos/rcasero/
--
Dr. Ramón Casero Cañas
Oxford e-Research Centre (OeRC)
University of Oxford
7 Keble Rd
Oxford OX1 3QG
tlf +44 (0) 1865 610739
web http://www.cs.ox.ac.uk/people/Ramon.CaseroCanas
photos http://www.flickr.com/photos/rcasero/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20131011/081bb842/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: mesh.png
Type: image/png
Size: 13847 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20131011/081bb842/attachment.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: slice231.png
Type: image/png
Size: 3524 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20131011/081bb842/attachment-0001.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: slice300.png
Type: image/png
Size: 3553 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20131011/081bb842/attachment-0002.png>
More information about the Insight-users
mailing list